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