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