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