1 : /******************************************************************************
2 : * $Id: gdal_crs.c 22670 2011-07-08 17:15:19Z rouault $
3 : *
4 : * Project: Mapinfo Image Warper
5 : * Purpose: Implemention of the GDALTransformer wrapper around CRS.C functions
6 : * to build a polynomial transformation based on ground control
7 : * points.
8 : * Author: Frank Warmerdam, warmerdam@pobox.com
9 : *
10 : ***************************************************************************
11 :
12 : CRS.C - Center for Remote Sensing rectification routines
13 :
14 : Written By: Brian J. Buckley
15 :
16 : At: The Center for Remote Sensing
17 : Michigan State University
18 : 302 Berkey Hall
19 : East Lansing, MI 48824
20 : (517)353-7195
21 :
22 : Written: 12/19/91
23 :
24 : Last Update: 12/26/91 Brian J. Buckley
25 : Last Update: 1/24/92 Brian J. Buckley
26 : Added printout of trnfile. Triggered by BDEBUG.
27 : Last Update: 1/27/92 Brian J. Buckley
28 : Fixed bug so that only the active control points were used.
29 : Last Update: 6/29/2011 C. F. Stallmann & R. van den Dool (South African National Space Agency)
30 : GCP refinement added
31 :
32 :
33 : Copyright (c) 1992, Michigan State University
34 :
35 : Permission is hereby granted, free of charge, to any person obtaining a
36 : copy of this software and associated documentation files (the "Software"),
37 : to deal in the Software without restriction, including without limitation
38 : the rights to use, copy, modify, merge, publish, distribute, sublicense,
39 : and/or sell copies of the Software, and to permit persons to whom the
40 : Software is furnished to do so, subject to the following conditions:
41 :
42 : The above copyright notice and this permission notice shall be included
43 : in all copies or substantial portions of the Software.
44 :
45 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
46 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
47 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
48 : THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
49 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
50 : FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
51 : DEALINGS IN THE SOFTWARE.
52 :
53 : ****************************************************************************/
54 :
55 : #include "gdal_alg.h"
56 : #include "cpl_conv.h"
57 : #include "cpl_minixml.h"
58 : #include "cpl_string.h"
59 :
60 : CPL_CVSID("$Id: gdal_crs.c 22670 2011-07-08 17:15:19Z rouault $");
61 :
62 : #define MAXORDER 3
63 :
64 : struct Control_Points
65 : {
66 : int count;
67 : double *e1;
68 : double *n1;
69 : double *e2;
70 : double *n2;
71 : int *status;
72 : };
73 :
74 : typedef struct
75 : {
76 : GDALTransformerInfo sTI;
77 :
78 : double adfToGeoX[20];
79 : double adfToGeoY[20];
80 :
81 : double adfFromGeoX[20];
82 : double adfFromGeoY[20];
83 :
84 : int nOrder;
85 : int bReversed;
86 :
87 : int nGCPCount;
88 : GDAL_GCP *pasGCPList;
89 : int bRefine;
90 : int nMinimumGcps;
91 : double dfTolerance;
92 :
93 : } GCPTransformInfo;
94 :
95 : CPL_C_START
96 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
97 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
98 : CPL_C_END
99 :
100 : /* crs.c */
101 : static int CRS_georef(double, double, double *, double *,
102 : double [], double [], int);
103 : static int CRS_compute_georef_equations(struct Control_Points *,
104 : double [], double [], double [], double [], int);
105 : static int remove_outliers(GCPTransformInfo *);
106 :
107 : static char *CRS_error_message[] = {
108 : "Failed to compute GCP transform: Not enough points available",
109 : "Failed to compute GCP transform: Transform is not solvable",
110 : "Failed to compute GCP transform: Not enough memory",
111 : "Failed to compute GCP transform: Parameter error",
112 : "Failed to compute GCP transform: Internal error"
113 : };
114 :
115 :
116 : /************************************************************************/
117 : /* GDALCreateGCPTransformer() */
118 : /************************************************************************/
119 :
120 : /**
121 : * Create GCP based polynomial transformer.
122 : *
123 : * Computes least squares fit polynomials from a provided set of GCPs,
124 : * and stores the coefficients for later transformation of points between
125 : * pixel/line and georeferenced coordinates.
126 : *
127 : * The return value should be used as a TransformArg in combination with
128 : * the transformation function GDALGCPTransform which fits the
129 : * GDALTransformerFunc signature. The returned transform argument should
130 : * be deallocated with GDALDestroyGCPTransformer when no longer needed.
131 : *
132 : * This function may fail (returning NULL) if the provided set of GCPs
133 : * are inadequate for the requested order, the determinate is zero or they
134 : * are otherwise "ill conditioned".
135 : *
136 : * Note that 2nd order requires at least 6 GCPs, and 3rd order requires at
137 : * least 10 gcps. If nReqOrder is 0 the highest order possible with the
138 : * provided gcp count will be used.
139 : *
140 : * @param nGCPCount the number of GCPs in pasGCPList.
141 : * @param pasGCPList an array of GCPs to be used as input.
142 : * @param nReqOrder the requested polynomial order. It should be 1, 2 or 3.
143 : *
144 : * @return the transform argument or NULL if creation fails.
145 : */
146 :
147 27 : void *GDALCreateGCPTransformerEx( int nGCPCount, const GDAL_GCP *pasGCPList,
148 : int nReqOrder, int bReversed, int bRefine, double dfTolerance, int nMinimumGcps)
149 :
150 : {
151 : GCPTransformInfo *psInfo;
152 : double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
153 : int *panStatus, iGCP;
154 : int nCRSresult;
155 : struct Control_Points sPoints;
156 :
157 27 : if( nReqOrder == 0 )
158 : {
159 20 : if( nGCPCount >= 10 )
160 0 : nReqOrder = 2; /*for now we avoid 3rd order since it is unstable*/
161 20 : else if( nGCPCount >= 6 )
162 1 : nReqOrder = 2;
163 : else
164 19 : nReqOrder = 1;
165 : }
166 :
167 27 : psInfo = (GCPTransformInfo *) CPLCalloc(sizeof(GCPTransformInfo),1);
168 27 : psInfo->bReversed = bReversed;
169 27 : psInfo->nOrder = nReqOrder;
170 27 : psInfo->bRefine = bRefine;
171 27 : psInfo->dfTolerance = dfTolerance;
172 27 : psInfo->nMinimumGcps = nMinimumGcps;
173 :
174 27 : psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
175 27 : psInfo->nGCPCount = nGCPCount;
176 :
177 27 : strcpy( psInfo->sTI.szSignature, "GTI" );
178 27 : psInfo->sTI.pszClassName = "GDALGCPTransformer";
179 27 : psInfo->sTI.pfnTransform = GDALGCPTransform;
180 27 : psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
181 27 : psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
182 :
183 : /* -------------------------------------------------------------------- */
184 : /* Compute the forward and reverse polynomials. */
185 : /* -------------------------------------------------------------------- */
186 :
187 27 : if(bRefine)
188 : {
189 1 : nCRSresult = remove_outliers(psInfo);
190 : }
191 : else
192 : {
193 : /* -------------------------------------------------------------------- */
194 : /* Allocate and initialize the working points list. */
195 : /* -------------------------------------------------------------------- */
196 26 : padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
197 26 : padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
198 26 : padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
199 26 : padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
200 26 : panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
201 :
202 130 : for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
203 : {
204 104 : panStatus[iGCP] = 1;
205 104 : padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
206 104 : padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
207 104 : padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
208 104 : padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
209 : }
210 :
211 26 : sPoints.count = nGCPCount;
212 26 : sPoints.e1 = padfRasterX;
213 26 : sPoints.n1 = padfRasterY;
214 26 : sPoints.e2 = padfGeoX;
215 26 : sPoints.n2 = padfGeoY;
216 26 : sPoints.status = panStatus;
217 26 : nCRSresult = CRS_compute_georef_equations( &sPoints,
218 : psInfo->adfToGeoX, psInfo->adfToGeoY,
219 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
220 : nReqOrder );
221 26 : CPLFree( padfGeoX );
222 26 : CPLFree( padfGeoY );
223 26 : CPLFree( padfRasterX );
224 26 : CPLFree( padfRasterY );
225 26 : CPLFree( panStatus );
226 : }
227 :
228 27 : if (nCRSresult != 1)
229 : {
230 0 : CPLError( CE_Failure, CPLE_AppDefined, "%s", CRS_error_message[-nCRSresult]);
231 0 : GDALDestroyGCPTransformer( psInfo );
232 0 : return NULL;
233 : }
234 : else
235 : {
236 27 : return psInfo;
237 : }
238 : }
239 :
240 26 : void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
241 : int nReqOrder, int bReversed )
242 :
243 : {
244 26 : return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, FALSE, -1, -1);
245 : }
246 :
247 1 : void *GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
248 : int nReqOrder, int bReversed, double dfTolerance, int nMinimumGcps)
249 :
250 : {
251 : //If no minimumGcp parameter was passed, we use the default value according to the model
252 1 : if(nMinimumGcps == -1)
253 : {
254 0 : nMinimumGcps = ((nReqOrder+1) * (nReqOrder+2)) / 2 + 1;
255 : }
256 1 : return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, dfTolerance, nMinimumGcps);
257 : }
258 :
259 :
260 : /************************************************************************/
261 : /* GDALDestroyGCPTransformer() */
262 : /************************************************************************/
263 :
264 : /**
265 : * Destroy GCP transformer.
266 : *
267 : * This function is used to destroy information about a GCP based
268 : * polynomial transformation created with GDALCreateGCPTransformer().
269 : *
270 : * @param pTransformArg the transform arg previously returned by
271 : * GDALCreateGCPTransformer().
272 : */
273 :
274 27 : void GDALDestroyGCPTransformer( void *pTransformArg )
275 :
276 : {
277 27 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
278 :
279 27 : VALIDATE_POINTER0( pTransformArg, "GDALDestroyGCPTransformer" );
280 :
281 27 : GDALDeinitGCPs( psInfo->nGCPCount, psInfo->pasGCPList );
282 27 : CPLFree( psInfo->pasGCPList );
283 :
284 27 : CPLFree( pTransformArg );
285 : }
286 :
287 : /************************************************************************/
288 : /* GDALGCPTransform() */
289 : /************************************************************************/
290 :
291 : /**
292 : * Transforms point based on GCP derived polynomial model.
293 : *
294 : * This function matches the GDALTransformerFunc signature, and can be
295 : * used to transform one or more points from pixel/line coordinates to
296 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
297 : *
298 : * @param pTransformArg return value from GDALCreateGCPTransformer().
299 : * @param bDstToSrc TRUE if transformation is from the destination
300 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
301 : * from pixel/line to georeferenced coordinates.
302 : * @param nPointCount the number of values in the x, y and z arrays.
303 : * @param x array containing the X values to be transformed.
304 : * @param y array containing the Y values to be transformed.
305 : * @param z array containing the Z values to be transformed.
306 : * @param panSuccess array in which a flag indicating success (TRUE) or
307 : * failure (FALSE) of the transformation are placed.
308 : *
309 : * @return TRUE.
310 : */
311 :
312 12223 : int GDALGCPTransform( void *pTransformArg, int bDstToSrc,
313 : int nPointCount,
314 : double *x, double *y, double *z,
315 : int *panSuccess )
316 :
317 : {
318 : int i;
319 12223 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
320 :
321 12223 : if( psInfo->bReversed )
322 0 : bDstToSrc = !bDstToSrc;
323 :
324 103106 : for( i = 0; i < nPointCount; i++ )
325 : {
326 90883 : if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
327 : {
328 0 : panSuccess[i] = FALSE;
329 0 : continue;
330 : }
331 :
332 90883 : if( bDstToSrc )
333 : {
334 81898 : CRS_georef( x[i], y[i], x + i, y + i,
335 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
336 : psInfo->nOrder );
337 : }
338 : else
339 : {
340 8985 : CRS_georef( x[i], y[i], x + i, y + i,
341 : psInfo->adfToGeoX, psInfo->adfToGeoY,
342 : psInfo->nOrder );
343 : }
344 90883 : panSuccess[i] = TRUE;
345 : }
346 :
347 12223 : return TRUE;
348 : }
349 :
350 : /************************************************************************/
351 : /* GDALSerializeGCPTransformer() */
352 : /************************************************************************/
353 :
354 3 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg )
355 :
356 : {
357 : CPLXMLNode *psTree;
358 3 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
359 :
360 3 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeGCPTransformer", NULL );
361 :
362 3 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "GCPTransformer" );
363 :
364 : /* -------------------------------------------------------------------- */
365 : /* Serialize Order and bReversed. */
366 : /* -------------------------------------------------------------------- */
367 3 : CPLCreateXMLElementAndValue(
368 : psTree, "Order",
369 : CPLSPrintf( "%d", psInfo->nOrder ) );
370 :
371 3 : CPLCreateXMLElementAndValue(
372 : psTree, "Reversed",
373 : CPLSPrintf( "%d", psInfo->bReversed ) );
374 :
375 3 : if( psInfo->bRefine )
376 : {
377 0 : CPLCreateXMLElementAndValue(
378 : psTree, "Refine",
379 : CPLSPrintf( "%d", psInfo->bRefine ) );
380 :
381 0 : CPLCreateXMLElementAndValue(
382 : psTree, "MinimumGcps",
383 : CPLSPrintf( "%d", psInfo->nMinimumGcps ) );
384 :
385 0 : CPLCreateXMLElementAndValue(
386 : psTree, "Tolerance",
387 : CPLSPrintf( "%f", psInfo->dfTolerance ) );
388 : }
389 :
390 : /* -------------------------------------------------------------------- */
391 : /* Attach GCP List. */
392 : /* -------------------------------------------------------------------- */
393 3 : if( psInfo->nGCPCount > 0 )
394 : {
395 : int iGCP;
396 3 : CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element,
397 3 : "GCPList" );
398 :
399 3 : if(psInfo->bRefine)
400 : {
401 0 : remove_outliers(psInfo);
402 : }
403 :
404 15 : for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
405 : {
406 : CPLXMLNode *psXMLGCP;
407 12 : GDAL_GCP *psGCP = psInfo->pasGCPList + iGCP;
408 :
409 12 : psXMLGCP = CPLCreateXMLNode( psGCPList, CXT_Element, "GCP" );
410 :
411 12 : CPLSetXMLValue( psXMLGCP, "#Id", psGCP->pszId );
412 :
413 12 : if( psGCP->pszInfo != NULL && strlen(psGCP->pszInfo) > 0 )
414 0 : CPLSetXMLValue( psXMLGCP, "Info", psGCP->pszInfo );
415 :
416 12 : CPLSetXMLValue( psXMLGCP, "#Pixel",
417 : CPLSPrintf( "%.4f", psGCP->dfGCPPixel ) );
418 :
419 12 : CPLSetXMLValue( psXMLGCP, "#Line",
420 : CPLSPrintf( "%.4f", psGCP->dfGCPLine ) );
421 :
422 12 : CPLSetXMLValue( psXMLGCP, "#X",
423 : CPLSPrintf( "%.12E", psGCP->dfGCPX ) );
424 :
425 12 : CPLSetXMLValue( psXMLGCP, "#Y",
426 : CPLSPrintf( "%.12E", psGCP->dfGCPY ) );
427 :
428 12 : if( psGCP->dfGCPZ != 0.0 )
429 0 : CPLSetXMLValue( psXMLGCP, "#GCPZ",
430 : CPLSPrintf( "%.12E", psGCP->dfGCPZ ) );
431 : }
432 : }
433 :
434 3 : return psTree;
435 : }
436 :
437 : /************************************************************************/
438 : /* GDALDeserializeReprojectionTransformer() */
439 : /************************************************************************/
440 :
441 5 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree )
442 :
443 : {
444 5 : GDAL_GCP *pasGCPList = 0;
445 5 : int nGCPCount = 0;
446 : void *pResult;
447 : int nReqOrder;
448 : int bReversed;
449 : int bRefine;
450 : int nMinimumGcps;
451 : double dfTolerance;
452 :
453 : /* -------------------------------------------------------------------- */
454 : /* Check for GCPs. */
455 : /* -------------------------------------------------------------------- */
456 5 : CPLXMLNode *psGCPList = CPLGetXMLNode( psTree, "GCPList" );
457 :
458 5 : if( psGCPList != NULL )
459 : {
460 5 : int nGCPMax = 0;
461 : CPLXMLNode *psXMLGCP;
462 :
463 : // Count GCPs.
464 36 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
465 26 : psXMLGCP = psXMLGCP->psNext )
466 26 : nGCPMax++;
467 :
468 5 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),nGCPMax);
469 :
470 36 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
471 26 : psXMLGCP = psXMLGCP->psNext )
472 : {
473 26 : GDAL_GCP *psGCP = pasGCPList + nGCPCount;
474 :
475 52 : if( !EQUAL(psXMLGCP->pszValue,"GCP") ||
476 26 : psXMLGCP->eType != CXT_Element )
477 0 : continue;
478 :
479 26 : GDALInitGCPs( 1, psGCP );
480 :
481 26 : CPLFree( psGCP->pszId );
482 26 : psGCP->pszId = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Id",""));
483 :
484 26 : CPLFree( psGCP->pszInfo );
485 26 : psGCP->pszInfo = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Info",""));
486 :
487 26 : psGCP->dfGCPPixel = atof(CPLGetXMLValue(psXMLGCP,"Pixel","0.0"));
488 26 : psGCP->dfGCPLine = atof(CPLGetXMLValue(psXMLGCP,"Line","0.0"));
489 :
490 26 : psGCP->dfGCPX = atof(CPLGetXMLValue(psXMLGCP,"X","0.0"));
491 26 : psGCP->dfGCPY = atof(CPLGetXMLValue(psXMLGCP,"Y","0.0"));
492 26 : psGCP->dfGCPZ = atof(CPLGetXMLValue(psXMLGCP,"Z","0.0"));
493 26 : nGCPCount++;
494 : }
495 : }
496 :
497 : /* -------------------------------------------------------------------- */
498 : /* Get other flags. */
499 : /* -------------------------------------------------------------------- */
500 5 : nReqOrder = atoi(CPLGetXMLValue(psTree,"Order","3"));
501 5 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
502 5 : bRefine = atoi(CPLGetXMLValue(psTree,"Refine","0"));
503 5 : nMinimumGcps = atoi(CPLGetXMLValue(psTree,"MinimumGcps","6"));
504 5 : dfTolerance = atof(CPLGetXMLValue(psTree,"Tolerance","1.0"));
505 :
506 : /* -------------------------------------------------------------------- */
507 : /* Generate transformation. */
508 : /* -------------------------------------------------------------------- */
509 5 : if(bRefine)
510 : {
511 1 : pResult = GDALCreateGCPRefineTransformer( nGCPCount, pasGCPList, nReqOrder,
512 : bReversed, dfTolerance, nMinimumGcps );
513 : }
514 : else
515 : {
516 4 : pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
517 : bReversed );
518 : }
519 :
520 : /* -------------------------------------------------------------------- */
521 : /* Cleanup GCP copy. */
522 : /* -------------------------------------------------------------------- */
523 5 : GDALDeinitGCPs( nGCPCount, pasGCPList );
524 5 : CPLFree( pasGCPList );
525 :
526 5 : return pResult;
527 : }
528 :
529 : /************************************************************************/
530 : /* ==================================================================== */
531 : /* Everything below this point derived from the CRS.C from GRASS. */
532 : /* ==================================================================== */
533 : /************************************************************************/
534 :
535 :
536 : /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT
537 : SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
538 :
539 : struct MATRIX
540 : {
541 : int n; /* SIZE OF THIS MATRIX (N x N) */
542 : double *v;
543 : };
544 :
545 : /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
546 :
547 : #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1]
548 :
549 :
550 : #define MSUCCESS 1 /* SUCCESS */
551 : #define MNPTERR 0 /* NOT ENOUGH POINTS */
552 : #define MUNSOLVABLE -1 /* NOT SOLVABLE */
553 : #define MMEMERR -2 /* NOT ENOUGH MEMORY */
554 : #define MPARMERR -3 /* PARAMETER ERROR */
555 : #define MINTERR -4 /* INTERNAL ERROR */
556 :
557 : /***************************************************************************/
558 : /*
559 : FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
560 : */
561 : /***************************************************************************/
562 :
563 : static int calccoef(struct Control_Points *,double *,double *,int);
564 : static int calcls(struct Control_Points *,struct MATRIX *,
565 : double *,double *,double *,double *);
566 : static int exactdet(struct Control_Points *,struct MATRIX *,
567 : double *,double *,double *,double *);
568 : static int solvemat(struct MATRIX *,double *,double *,double *,double *);
569 : static double term(int,double,double);
570 :
571 : /***************************************************************************/
572 : /*
573 : TRANSFORM A SINGLE COORDINATE PAIR.
574 : */
575 : /***************************************************************************/
576 :
577 : static int
578 90883 : CRS_georef (
579 : double e1, /* EASTINGS TO BE TRANSFORMED */
580 : double n1, /* NORTHINGS TO BE TRANSFORMED */
581 : double *e, /* EASTINGS TO BE TRANSFORMED */
582 : double *n, /* NORTHINGS TO BE TRANSFORMED */
583 : double E[], /* EASTING COEFFICIENTS */
584 : double N[], /* NORTHING COEFFICIENTS */
585 : int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
586 : ORDER USED TO CALCULATE THE COEFFICIENTS */
587 : )
588 : {
589 : double e3, e2n, en2, n3, e2, en, n2;
590 :
591 90883 : switch(order)
592 : {
593 : case 1:
594 :
595 88739 : *e = E[0] + E[1] * e1 + E[2] * n1;
596 88739 : *n = N[0] + N[1] * e1 + N[2] * n1;
597 88739 : break;
598 :
599 : case 2:
600 :
601 2144 : e2 = e1 * e1;
602 2144 : n2 = n1 * n1;
603 2144 : en = e1 * n1;
604 :
605 4288 : *e = E[0] + E[1] * e1 + E[2] * n1 +
606 2144 : E[3] * e2 + E[4] * en + E[5] * n2;
607 4288 : *n = N[0] + N[1] * e1 + N[2] * n1 +
608 2144 : N[3] * e2 + N[4] * en + N[5] * n2;
609 2144 : break;
610 :
611 : case 3:
612 :
613 0 : e2 = e1 * e1;
614 0 : en = e1 * n1;
615 0 : n2 = n1 * n1;
616 0 : e3 = e1 * e2;
617 0 : e2n = e2 * n1;
618 0 : en2 = e1 * n2;
619 0 : n3 = n1 * n2;
620 :
621 0 : *e = E[0] +
622 0 : E[1] * e1 + E[2] * n1 +
623 0 : E[3] * e2 + E[4] * en + E[5] * n2 +
624 0 : E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
625 0 : *n = N[0] +
626 0 : N[1] * e1 + N[2] * n1 +
627 0 : N[3] * e2 + N[4] * en + N[5] * n2 +
628 0 : N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
629 0 : break;
630 :
631 : default:
632 :
633 0 : return(MPARMERR);
634 : }
635 :
636 90883 : return(MSUCCESS);
637 : }
638 :
639 : /***************************************************************************/
640 : /*
641 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
642 : */
643 : /***************************************************************************/
644 :
645 : static int
646 29 : CRS_compute_georef_equations (struct Control_Points *cp,
647 : double E12[], double N12[],
648 : double E21[], double N21[],
649 : int order)
650 : {
651 : double *tempptr;
652 : int status;
653 :
654 29 : if(order < 1 || order > MAXORDER)
655 0 : return(MPARMERR);
656 :
657 : /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
658 :
659 29 : status = calccoef(cp,E12,N12,order);
660 29 : if(status != MSUCCESS)
661 0 : return(status);
662 :
663 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
664 :
665 29 : tempptr = cp->e1;
666 29 : cp->e1 = cp->e2;
667 29 : cp->e2 = tempptr;
668 29 : tempptr = cp->n1;
669 29 : cp->n1 = cp->n2;
670 29 : cp->n2 = tempptr;
671 :
672 : /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
673 :
674 29 : status = calccoef(cp,E21,N21,order);
675 :
676 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
677 :
678 29 : tempptr = cp->e1;
679 29 : cp->e1 = cp->e2;
680 29 : cp->e2 = tempptr;
681 29 : tempptr = cp->n1;
682 29 : cp->n1 = cp->n2;
683 29 : cp->n2 = tempptr;
684 :
685 29 : return(status);
686 : }
687 :
688 : /***************************************************************************/
689 : /*
690 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
691 : */
692 : /***************************************************************************/
693 :
694 : static int
695 58 : calccoef (struct Control_Points *cp, double E[], double N[], int order)
696 : {
697 : struct MATRIX m;
698 : double *a;
699 : double *b;
700 : int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
701 : int status, i;
702 :
703 : /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
704 :
705 320 : for(i = numactive = 0 ; i < cp->count ; i++)
706 : {
707 262 : if(cp->status[i] > 0)
708 262 : numactive++;
709 : }
710 :
711 : /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
712 : A TRANSFORMATION OF THIS ORDER */
713 :
714 58 : m.n = ((order + 1) * (order + 2)) / 2;
715 :
716 58 : if(numactive < m.n)
717 0 : return(MNPTERR);
718 :
719 : /* INITIALIZE MATRIX */
720 :
721 58 : m.v = (double *)CPLCalloc(m.n*m.n,sizeof(double));
722 58 : if(m.v == NULL)
723 : {
724 0 : return(MMEMERR);
725 : }
726 58 : a = (double *)CPLCalloc(m.n,sizeof(double));
727 58 : if(a == NULL)
728 : {
729 0 : CPLFree((char *)m.v);
730 0 : return(MMEMERR);
731 : }
732 58 : b = (double *)CPLCalloc(m.n,sizeof(double));
733 58 : if(b == NULL)
734 : {
735 0 : CPLFree((char *)m.v);
736 0 : CPLFree((char *)a);
737 0 : return(MMEMERR);
738 : }
739 :
740 58 : if(numactive == m.n)
741 8 : status = exactdet(cp,&m,a,b,E,N);
742 : else
743 50 : status = calcls(cp,&m,a,b,E,N);
744 :
745 58 : CPLFree((char *)m.v);
746 58 : CPLFree((char *)a);
747 58 : CPLFree((char *)b);
748 :
749 58 : return(status);
750 : }
751 :
752 : /***************************************************************************/
753 : /*
754 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
755 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
756 : */
757 : /***************************************************************************/
758 :
759 8 : static int exactdet (
760 : struct Control_Points *cp,
761 : struct MATRIX *m,
762 : double a[],
763 : double b[],
764 : double E[], /* EASTING COEFFICIENTS */
765 : double N[] /* NORTHING COEFFICIENTS */
766 : )
767 : {
768 : int pntnow, currow, j;
769 :
770 8 : currow = 1;
771 32 : for(pntnow = 0 ; pntnow < cp->count ; pntnow++)
772 : {
773 24 : if(cp->status[pntnow] > 0)
774 : {
775 : /* POPULATE MATRIX M */
776 :
777 96 : for(j = 1 ; j <= m->n ; j++)
778 : {
779 72 : M(currow,j) = term(j,cp->e1[pntnow],cp->n1[pntnow]);
780 : }
781 :
782 : /* POPULATE MATRIX A AND B */
783 :
784 24 : a[currow-1] = cp->e2[pntnow];
785 24 : b[currow-1] = cp->n2[pntnow];
786 :
787 24 : currow++;
788 : }
789 : }
790 :
791 8 : if(currow - 1 != m->n)
792 0 : return(MINTERR);
793 :
794 8 : return(solvemat(m,a,b,E,N));
795 : }
796 :
797 : /***************************************************************************/
798 : /*
799 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
800 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS
801 : ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
802 : */
803 : /***************************************************************************/
804 :
805 50 : static int calcls (
806 : struct Control_Points *cp,
807 : struct MATRIX *m,
808 : double a[],
809 : double b[],
810 : double E[], /* EASTING COEFFICIENTS */
811 : double N[] /* NORTHING COEFFICIENTS */
812 : )
813 : {
814 50 : int i, j, n, numactive = 0;
815 :
816 : /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
817 :
818 206 : for(i = 1 ; i <= m->n ; i++)
819 : {
820 486 : for(j = i ; j <= m->n ; j++)
821 330 : M(i,j) = 0.0;
822 156 : a[i-1] = b[i-1] = 0.0;
823 : }
824 :
825 : /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
826 : THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
827 :
828 288 : for(n = 0 ; n < cp->count ; n++)
829 : {
830 238 : if(cp->status[n] > 0)
831 : {
832 238 : numactive++;
833 1000 : for(i = 1 ; i <= m->n ; i++)
834 : {
835 2430 : for(j = i ; j <= m->n ; j++)
836 1668 : M(i,j) += term(i,cp->e1[n],cp->n1[n]) * term(j,cp->e1[n],cp->n1[n]);
837 :
838 762 : a[i-1] += cp->e2[n] * term(i,cp->e1[n],cp->n1[n]);
839 762 : b[i-1] += cp->n2[n] * term(i,cp->e1[n],cp->n1[n]);
840 : }
841 : }
842 : }
843 :
844 50 : if(numactive <= m->n)
845 0 : return(MINTERR);
846 :
847 : /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
848 :
849 156 : for(i = 2 ; i <= m->n ; i++)
850 : {
851 280 : for(j = 1 ; j < i ; j++)
852 174 : M(i,j) = M(j,i);
853 : }
854 :
855 50 : return(solvemat(m,a,b,E,N));
856 : }
857 :
858 : /***************************************************************************/
859 : /*
860 : CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
861 :
862 : ORDER\TERM 1 2 3 4 5 6 7 8 9 10
863 : 1 e0n0 e1n0 e0n1
864 : 2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
865 : 3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
866 : */
867 : /***************************************************************************/
868 :
869 4932 : static double term (int term, double e, double n)
870 : {
871 4932 : switch(term)
872 : {
873 1500 : case 1: return((double)1.0);
874 1500 : case 2: return((double)e);
875 1500 : case 3: return((double)n);
876 144 : case 4: return((double)(e*e));
877 144 : case 5: return((double)(e*n));
878 144 : case 6: return((double)(n*n));
879 0 : case 7: return((double)(e*e*e));
880 0 : case 8: return((double)(e*e*n));
881 0 : case 9: return((double)(e*n*n));
882 0 : case 10: return((double)(n*n*n));
883 : }
884 0 : return((double)0.0);
885 : }
886 :
887 : /***************************************************************************/
888 : /*
889 : SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
890 : GAUSSIAN ELIMINATION METHOD.
891 :
892 : | M11 M12 ... M1n | | E0 | | a0 |
893 : | M21 M22 ... M2n | | E1 | = | a1 |
894 : | . . . . | | . | | . |
895 : | Mn1 Mn2 ... Mnn | | En-1 | | an-1 |
896 :
897 : and
898 :
899 : | M11 M12 ... M1n | | N0 | | b0 |
900 : | M21 M22 ... M2n | | N1 | = | b1 |
901 : | . . . . | | . | | . |
902 : | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 |
903 : */
904 : /***************************************************************************/
905 :
906 58 : static int solvemat (struct MATRIX *m,
907 : double a[], double b[], double E[], double N[])
908 : {
909 : int i, j, i2, j2, imark;
910 : double factor, temp;
911 : double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
912 :
913 238 : for(i = 1 ; i <= m->n ; i++)
914 : {
915 180 : j = i;
916 :
917 : /* find row with largest magnitude value for pivot value */
918 :
919 180 : pivot = M(i,j);
920 180 : imark = i;
921 378 : for(i2 = i + 1 ; i2 <= m->n ; i2++)
922 : {
923 198 : temp = fabs(M(i2,j));
924 198 : if(temp > fabs(pivot))
925 : {
926 120 : pivot = M(i2,j);
927 120 : imark = i2;
928 : }
929 : }
930 :
931 : /* if the pivot is very small then the points are nearly co-linear */
932 : /* co-linear points result in an undefined matrix, and nearly */
933 : /* co-linear points results in a solution with rounding error */
934 :
935 180 : if(pivot == 0.0)
936 0 : return(MUNSOLVABLE);
937 :
938 : /* if row with highest pivot is not the current row, switch them */
939 :
940 180 : if(imark != i)
941 : {
942 396 : for(j2 = 1 ; j2 <= m->n ; j2++)
943 : {
944 303 : temp = M(imark,j2);
945 303 : M(imark,j2) = M(i,j2);
946 303 : M(i,j2) = temp;
947 : }
948 :
949 93 : temp = a[imark-1];
950 93 : a[imark-1] = a[i-1];
951 93 : a[i-1] = temp;
952 :
953 93 : temp = b[imark-1];
954 93 : b[imark-1] = b[i-1];
955 93 : b[i-1] = temp;
956 : }
957 :
958 : /* compute zeros above and below the pivot, and compute
959 : values for the rest of the row as well */
960 :
961 756 : for(i2 = 1 ; i2 <= m->n ; i2++)
962 : {
963 576 : if(i2 != i)
964 : {
965 396 : factor = M(i2,j) / pivot;
966 1278 : for(j2 = j ; j2 <= m->n ; j2++)
967 882 : M(i2,j2) -= factor * M(i,j2);
968 396 : a[i2-1] -= factor * a[i-1];
969 396 : b[i2-1] -= factor * b[i-1];
970 : }
971 : }
972 : }
973 :
974 : /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
975 : COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
976 :
977 238 : for(i = 1 ; i <= m->n ; i++)
978 : {
979 180 : E[i-1] = a[i-1] / M(i,i);
980 180 : N[i-1] = b[i-1] / M(i,i);
981 : }
982 :
983 58 : return(MSUCCESS);
984 : }
985 :
986 : /***************************************************************************/
987 : /*
988 : DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
989 : OUTLIER.
990 :
991 : THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
992 : AND THE ALLOWED TOLERANCE:
993 :
994 : sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
995 : lineAdj = b0 + b1*sample + b2*line + b3*line*sample
996 :
997 : WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
998 :
999 : [residualSample] = [A1][sampleCoefficients] - [b1]
1000 : [residualLine] = [A2][lineCoefficients] - [b2]
1001 :
1002 : sampleResidual^2 = sum( [residualSample]^2 )
1003 : lineResidual^2 = sum( [lineSample]^2 )
1004 :
1005 : residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
1006 :
1007 : THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
1008 : CONSIDERED THE WORST OUTLIER.
1009 :
1010 : IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
1011 : */
1012 : /***************************************************************************/
1013 3 : static int worst_outlier(struct Control_Points *cp, double E[], double N[], double dfTolerance)
1014 : {
1015 : double *padfResiduals;
1016 : int nI, nIndex;
1017 : double dfDifference, dfSampleResidual, dfLineResidual, dfSampleRes, dfLineRes, dfCurrentDifference;
1018 : double dfE1, dfN1, dfE2, dfN2, dfEn;
1019 :
1020 3 : padfResiduals = (double *) CPLCalloc(sizeof(double),cp->count);
1021 3 : dfSampleResidual = 0.0;
1022 3 : dfLineResidual = 0.0;
1023 :
1024 30 : for(nI = 0; nI < cp->count; nI++)
1025 : {
1026 27 : dfE1 = cp->e1[nI];
1027 27 : dfN1 = cp->n1[nI];
1028 27 : dfE2 = dfE1 * dfE1;
1029 27 : dfN2 = dfN1 * dfN1;
1030 27 : dfEn = dfE1 * dfN1;
1031 :
1032 27 : dfSampleRes = E[0] + E[1] * dfE1 + E[2] * dfN1 + E[3] * dfE2 + E[4] * dfEn + E[5] * dfN2 - cp->e2[nI];
1033 27 : dfLineRes = N[0] + N[1] * dfE1 + N[2] * dfN1 + N[3] * dfE2 + N[4] * dfEn + N[5] * dfN2 - cp->n2[nI];
1034 :
1035 27 : dfSampleResidual += dfSampleRes*dfSampleRes;
1036 27 : dfLineResidual += dfLineRes*dfLineRes;
1037 :
1038 27 : padfResiduals[nI] = sqrt(dfSampleRes*dfSampleRes + dfLineRes*dfLineRes);
1039 : }
1040 :
1041 3 : nIndex = -1;
1042 3 : dfDifference = -1.0;
1043 30 : for(nI = 0; nI < cp->count; nI++)
1044 : {
1045 27 : dfCurrentDifference = padfResiduals[nI];
1046 27 : if(fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
1047 : {
1048 8 : dfCurrentDifference = 0.0;
1049 : }
1050 27 : if(dfCurrentDifference > dfDifference && dfCurrentDifference >= dfTolerance)
1051 : {
1052 6 : dfDifference = dfCurrentDifference;
1053 6 : nIndex = nI;
1054 : }
1055 : }
1056 3 : CPLFree( padfResiduals );
1057 3 : return nIndex;
1058 : }
1059 :
1060 : /***************************************************************************/
1061 : /*
1062 : REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
1063 : ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
1064 :
1065 : 1. WE CALCULATE THE COEFFICIENTS FOR ALL THE GCPS.
1066 : 2. THE GCP LIST WILL BE SCANED TO DETERMINE THE WORST OUTLIER USING
1067 : THE CALCULATED COEFFICIENTS.
1068 : 3. THE WORST OUTLIER WILL BE REMOVED FROM THE GCP LIST.
1069 : 4. THE COEFFICIENTS WILL BE RECALCULATED WITHOUT THE WORST OUTLIER.
1070 : 5. STEP 1 TO 4 ARE EXECUTED UNTIL THE MINIMUM NUMBER OF GCPS WERE REACHED
1071 : OR IF NO GCP IS CONSIDERED AN OUTLIER WITH THE PASSED TOLERANCE.
1072 : */
1073 : /***************************************************************************/
1074 1 : static int remove_outliers( GCPTransformInfo *psInfo )
1075 : {
1076 : double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
1077 : int *panStatus;
1078 : int nI, nCRSresult, nGCPCount, nMinimumGcps, nReqOrder;
1079 : double dfTolerance;
1080 : struct Control_Points sPoints;
1081 :
1082 1 : nGCPCount = psInfo->nGCPCount;
1083 1 : nMinimumGcps = psInfo->nMinimumGcps;
1084 1 : nReqOrder = psInfo->nOrder;
1085 1 : dfTolerance = psInfo->dfTolerance;
1086 :
1087 1 : padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
1088 1 : padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
1089 1 : padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
1090 1 : padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
1091 1 : panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
1092 :
1093 11 : for( nI = 0; nI < nGCPCount; nI++ )
1094 : {
1095 10 : panStatus[nI] = 1;
1096 10 : padfGeoX[nI] = psInfo->pasGCPList[nI].dfGCPX;
1097 10 : padfGeoY[nI] = psInfo->pasGCPList[nI].dfGCPY;
1098 10 : padfRasterX[nI] = psInfo->pasGCPList[nI].dfGCPPixel;
1099 10 : padfRasterY[nI] = psInfo->pasGCPList[nI].dfGCPLine;
1100 : }
1101 :
1102 1 : sPoints.count = nGCPCount;
1103 1 : sPoints.e1 = padfRasterX;
1104 1 : sPoints.n1 = padfRasterY;
1105 1 : sPoints.e2 = padfGeoX;
1106 1 : sPoints.n2 = padfGeoY;
1107 1 : sPoints.status = panStatus;
1108 :
1109 1 : nCRSresult = CRS_compute_georef_equations( &sPoints,
1110 : psInfo->adfToGeoX, psInfo->adfToGeoY,
1111 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
1112 : nReqOrder );
1113 :
1114 4 : while(sPoints.count > nMinimumGcps)
1115 : {
1116 : int nIndex;
1117 :
1118 3 : nIndex = worst_outlier(&sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, dfTolerance);
1119 :
1120 : //If no outliers were detected, stop the GCP elimination
1121 3 : if(nIndex == -1)
1122 : {
1123 1 : break;
1124 : }
1125 :
1126 2 : CPLFree(psInfo->pasGCPList[nIndex].pszId);
1127 2 : CPLFree(psInfo->pasGCPList[nIndex].pszInfo);
1128 :
1129 8 : for( nI = nIndex; nI < sPoints.count - 1; nI++ )
1130 : {
1131 6 : sPoints.e1[nI] = sPoints.e1[nI + 1];
1132 6 : sPoints.n1[nI] = sPoints.n1[nI + 1];
1133 6 : sPoints.e2[nI] = sPoints.e2[nI + 1];
1134 6 : sPoints.n2[nI] = sPoints.n2[nI + 1];
1135 6 : psInfo->pasGCPList[nI].pszId = psInfo->pasGCPList[nI + 1].pszId;
1136 6 : psInfo->pasGCPList[nI].pszInfo = psInfo->pasGCPList[nI + 1].pszInfo;
1137 : }
1138 :
1139 2 : sPoints.count = sPoints.count - 1;
1140 :
1141 2 : nCRSresult = CRS_compute_georef_equations( &sPoints,
1142 : psInfo->adfToGeoX, psInfo->adfToGeoY,
1143 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
1144 : nReqOrder );
1145 : }
1146 :
1147 9 : for( nI = 0; nI < sPoints.count; nI++ )
1148 : {
1149 8 : psInfo->pasGCPList[nI].dfGCPX = sPoints.e2[nI];
1150 8 : psInfo->pasGCPList[nI].dfGCPY = sPoints.n2[nI];
1151 8 : psInfo->pasGCPList[nI].dfGCPPixel = sPoints.e1[nI];
1152 8 : psInfo->pasGCPList[nI].dfGCPLine = sPoints.n1[nI];
1153 : }
1154 1 : psInfo->nGCPCount = sPoints.count;
1155 :
1156 1 : CPLFree( sPoints.e1 );
1157 1 : CPLFree( sPoints.n1 );
1158 1 : CPLFree( sPoints.e2 );
1159 1 : CPLFree( sPoints.n2 );
1160 1 : CPLFree( sPoints.status );
1161 1 : return nCRSresult;
1162 : }
|
