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 if all points have been successfully transformed.
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 17812 : int bRet = TRUE;
440 58807 : for (i = 0; i < nPointCount; i++)
441 : {
442 40995 : 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 40995 : if (bDstToSrc)
450 : {
451 28341 : CRS_georef(x[i] - psInfo->x2_mean, y[i] - psInfo->y2_mean, x + i,
452 28341 : y + i, psInfo->adfFromGeoX, psInfo->adfFromGeoY,
453 : psInfo->nOrder);
454 : }
455 : else
456 : {
457 12654 : CRS_georef(x[i] - psInfo->x1_mean, y[i] - psInfo->y1_mean, x + i,
458 12654 : y + i, psInfo->adfToGeoX, psInfo->adfToGeoY,
459 : psInfo->nOrder);
460 : }
461 40995 : panSuccess[i] = TRUE;
462 : }
463 :
464 17812 : return bRet;
465 : }
466 :
467 : /************************************************************************/
468 : /* GDALSerializeGCPTransformer() */
469 : /************************************************************************/
470 :
471 9 : CPLXMLNode *GDALSerializeGCPTransformer(void *pTransformArg)
472 :
473 : {
474 9 : CPLXMLNode *psTree = nullptr;
475 9 : GCPTransformInfo *psInfo = static_cast<GCPTransformInfo *>(pTransformArg);
476 :
477 9 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeGCPTransformer", nullptr);
478 :
479 9 : psTree = CPLCreateXMLNode(nullptr, CXT_Element, "GCPTransformer");
480 :
481 : /* -------------------------------------------------------------------- */
482 : /* Serialize Order and bReversed. */
483 : /* -------------------------------------------------------------------- */
484 9 : CPLCreateXMLElementAndValue(psTree, "Order",
485 : CPLSPrintf("%d", psInfo->nOrder));
486 :
487 9 : CPLCreateXMLElementAndValue(psTree, "Reversed",
488 : CPLSPrintf("%d", psInfo->bReversed));
489 :
490 9 : 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 9 : if (!psInfo->asGCPs.empty())
506 : {
507 9 : if (psInfo->bRefine)
508 : {
509 0 : remove_outliers(psInfo);
510 : }
511 :
512 9 : GDALSerializeGCPListToXML(psTree, psInfo->asGCPs, nullptr);
513 : }
514 :
515 9 : return psTree;
516 : }
517 :
518 : /************************************************************************/
519 : /* GDALDeserializeReprojectionTransformer() */
520 : /************************************************************************/
521 :
522 12 : void *GDALDeserializeGCPTransformer(CPLXMLNode *psTree)
523 :
524 : {
525 12 : std::vector<gdal::GCP> asGCPs;
526 12 : void *pResult = nullptr;
527 12 : int nReqOrder = 0;
528 12 : int bReversed = 0;
529 12 : int bRefine = 0;
530 12 : int nMinimumGcps = 0;
531 12 : double dfTolerance = 0.0;
532 :
533 : /* -------------------------------------------------------------------- */
534 : /* Check for GCPs. */
535 : /* -------------------------------------------------------------------- */
536 12 : CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList");
537 :
538 12 : if (psGCPList != nullptr)
539 : {
540 12 : GDALDeserializeGCPListFromXML(psGCPList, asGCPs, nullptr);
541 : }
542 :
543 : /* -------------------------------------------------------------------- */
544 : /* Get other flags. */
545 : /* -------------------------------------------------------------------- */
546 12 : nReqOrder = atoi(CPLGetXMLValue(psTree, "Order", "3"));
547 12 : bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
548 12 : bRefine = atoi(CPLGetXMLValue(psTree, "Refine", "0"));
549 12 : nMinimumGcps = atoi(CPLGetXMLValue(psTree, "MinimumGcps", "6"));
550 12 : dfTolerance = CPLAtof(CPLGetXMLValue(psTree, "Tolerance", "1.0"));
551 :
552 : /* -------------------------------------------------------------------- */
553 : /* Generate transformation. */
554 : /* -------------------------------------------------------------------- */
555 12 : 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 11 : pResult = GDALCreateGCPTransformer(static_cast<int>(asGCPs.size()),
564 : gdal::GCP::c_ptr(asGCPs), nReqOrder,
565 : bReversed);
566 : }
567 :
568 24 : 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 : /*
592 : FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
593 : */
594 : /***************************************************************************/
595 :
596 : static int calccoef(struct Control_Points *, double, double, double *, double *,
597 : int);
598 : static int calcls(struct Control_Points *, struct MATRIX *, double, double,
599 : double *, double *, double *, double *);
600 : static int exactdet(struct Control_Points *, struct MATRIX *, double, double,
601 : double *, double *, double *, double *);
602 : static int solvemat(struct MATRIX *, double *, double *, double *, double *);
603 : static double term(int, double, double);
604 :
605 : /***************************************************************************/
606 : /*
607 : TRANSFORM A SINGLE COORDINATE PAIR.
608 : */
609 : /***************************************************************************/
610 :
611 : static int
612 41022 : CRS_georef(double e1, /* EASTINGS TO BE TRANSFORMED */
613 : double n1, /* NORTHINGS TO BE TRANSFORMED */
614 : double *e, /* EASTINGS TO BE TRANSFORMED */
615 : double *n, /* NORTHINGS TO BE TRANSFORMED */
616 : double E[], /* EASTING COEFFICIENTS */
617 : double N[], /* NORTHING COEFFICIENTS */
618 : int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
619 : ORDER USED TO CALCULATE THE COEFFICIENTS */
620 : )
621 : {
622 41022 : double e3 = 0.0;
623 41022 : double e2n = 0.0;
624 41022 : double en2 = 0.0;
625 41022 : double n3 = 0.0;
626 41022 : double e2 = 0.0;
627 41022 : double en = 0.0;
628 41022 : double n2 = 0.0;
629 :
630 41022 : switch (order)
631 : {
632 38869 : case 1:
633 :
634 38869 : *e = E[0] + E[1] * e1 + E[2] * n1;
635 38869 : *n = N[0] + N[1] * e1 + N[2] * n1;
636 38869 : break;
637 :
638 2153 : case 2:
639 :
640 2153 : e2 = e1 * e1;
641 2153 : n2 = n1 * n1;
642 2153 : en = e1 * n1;
643 :
644 2153 : *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en +
645 2153 : E[5] * n2;
646 2153 : *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en +
647 2153 : N[5] * n2;
648 2153 : break;
649 :
650 0 : case 3:
651 :
652 0 : e2 = e1 * e1;
653 0 : en = e1 * n1;
654 0 : n2 = n1 * n1;
655 0 : e3 = e1 * e2;
656 0 : e2n = e2 * n1;
657 0 : en2 = e1 * n2;
658 0 : n3 = n1 * n2;
659 :
660 0 : *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en +
661 0 : E[5] * n2 + E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
662 0 : *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en +
663 0 : N[5] * n2 + N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
664 0 : break;
665 :
666 0 : default:
667 :
668 0 : return (MPARMERR);
669 : }
670 :
671 41022 : return (MSUCCESS);
672 : }
673 :
674 : /***************************************************************************/
675 : /*
676 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
677 : */
678 : /***************************************************************************/
679 :
680 53 : static int CRS_compute_georef_equations(GCPTransformInfo *psInfo,
681 : struct Control_Points *cp, double E12[],
682 : double N12[], double E21[],
683 : double N21[], int order)
684 : {
685 53 : double *tempptr = nullptr;
686 53 : int status = 0;
687 :
688 53 : if (order < 1 || order > MAXORDER)
689 0 : return (MPARMERR);
690 :
691 : /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
692 :
693 53 : status = calccoef(cp, psInfo->x1_mean, psInfo->y1_mean, E12, N12, order);
694 53 : if (status != MSUCCESS)
695 1 : return (status);
696 :
697 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
698 :
699 52 : tempptr = cp->e1;
700 52 : cp->e1 = cp->e2;
701 52 : cp->e2 = tempptr;
702 52 : tempptr = cp->n1;
703 52 : cp->n1 = cp->n2;
704 52 : cp->n2 = tempptr;
705 :
706 : /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
707 :
708 52 : status = calccoef(cp, psInfo->x2_mean, psInfo->y2_mean, E21, N21, order);
709 :
710 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
711 :
712 52 : tempptr = cp->e1;
713 52 : cp->e1 = cp->e2;
714 52 : cp->e2 = tempptr;
715 52 : tempptr = cp->n1;
716 52 : cp->n1 = cp->n2;
717 52 : cp->n2 = tempptr;
718 :
719 52 : return (status);
720 : }
721 :
722 : /***************************************************************************/
723 : /*
724 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
725 : */
726 : /***************************************************************************/
727 :
728 105 : static int calccoef(struct Control_Points *cp, double x_mean, double y_mean,
729 : double E[], double N[], int order)
730 : {
731 : struct MATRIX m;
732 105 : double *a = nullptr;
733 105 : double *b = nullptr;
734 105 : int numactive = 0; /* NUMBER OF ACTIVE CONTROL POINTS */
735 105 : int status = 0;
736 105 : int i = 0;
737 :
738 105 : memset(&m, 0, sizeof(m));
739 :
740 : /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
741 :
742 1772 : for (i = numactive = 0; i < cp->count; i++)
743 : {
744 1667 : if (cp->status[i] > 0)
745 1667 : numactive++;
746 : }
747 :
748 : /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
749 : A TRANSFORMATION OF THIS ORDER */
750 :
751 105 : m.n = ((order + 1) * (order + 2)) / 2;
752 :
753 105 : if (numactive < m.n)
754 1 : return (MNPTERR);
755 :
756 : /* INITIALIZE MATRIX */
757 :
758 104 : m.v = static_cast<double *>(
759 104 : VSICalloc(cpl::fits_on<int>(m.n * m.n), sizeof(double)));
760 104 : if (m.v == nullptr)
761 : {
762 0 : return (MMEMERR);
763 : }
764 104 : a = static_cast<double *>(VSICalloc(m.n, sizeof(double)));
765 104 : if (a == nullptr)
766 : {
767 0 : CPLFree(m.v);
768 0 : return (MMEMERR);
769 : }
770 104 : b = static_cast<double *>(VSICalloc(m.n, sizeof(double)));
771 104 : if (b == nullptr)
772 : {
773 0 : CPLFree(m.v);
774 0 : CPLFree(a);
775 0 : return (MMEMERR);
776 : }
777 :
778 104 : if (numactive == m.n)
779 24 : status = exactdet(cp, &m, x_mean, y_mean, a, b, E, N);
780 : else
781 80 : status = calcls(cp, &m, x_mean, y_mean, a, b, E, N);
782 :
783 104 : CPLFree(m.v);
784 104 : CPLFree(a);
785 104 : CPLFree(b);
786 :
787 104 : return (status);
788 : }
789 :
790 : /***************************************************************************/
791 : /*
792 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
793 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
794 : */
795 : /***************************************************************************/
796 :
797 24 : static int exactdet(struct Control_Points *cp, struct MATRIX *m, double x_mean,
798 : double y_mean, double a[], double b[],
799 : double E[], /* EASTING COEFFICIENTS */
800 : double N[] /* NORTHING COEFFICIENTS */
801 : )
802 : {
803 24 : int currow = 1;
804 :
805 96 : for (int pntnow = 0; pntnow < cp->count; pntnow++)
806 : {
807 72 : if (cp->status[pntnow] > 0)
808 : {
809 : /* POPULATE MATRIX M */
810 :
811 288 : for (int j = 1; j <= m->n; j++)
812 : {
813 216 : M(currow, j) =
814 216 : term(j, cp->e1[pntnow] - x_mean, cp->n1[pntnow] - y_mean);
815 : }
816 :
817 : /* POPULATE MATRIX A AND B */
818 :
819 72 : a[currow - 1] = cp->e2[pntnow];
820 72 : b[currow - 1] = cp->n2[pntnow];
821 :
822 72 : currow++;
823 : }
824 : }
825 :
826 24 : if (currow - 1 != m->n)
827 0 : return (MINTERR);
828 :
829 24 : return (solvemat(m, a, b, E, N));
830 : }
831 :
832 : /***************************************************************************/
833 : /*
834 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
835 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS
836 : ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
837 : */
838 : /***************************************************************************/
839 :
840 80 : static int calcls(struct Control_Points *cp, struct MATRIX *m, double x_mean,
841 : double y_mean, double a[], double b[],
842 : double E[], /* EASTING COEFFICIENTS */
843 : double N[] /* NORTHING COEFFICIENTS */
844 : )
845 : {
846 80 : int numactive = 0;
847 :
848 : /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
849 :
850 344 : for (int i = 1; i <= m->n; i++)
851 : {
852 864 : for (int j = i; j <= m->n; j++)
853 600 : M(i, j) = 0.0;
854 264 : a[i - 1] = b[i - 1] = 0.0;
855 : }
856 :
857 : /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
858 : THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
859 :
860 1674 : for (int n = 0; n < cp->count; n++)
861 : {
862 1594 : if (cp->status[n] > 0)
863 : {
864 1594 : numactive++;
865 10204 : for (int i = 1; i <= m->n; i++)
866 : {
867 37314 : for (int j = i; j <= m->n; j++)
868 57408 : M(i, j) += term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean) *
869 28704 : term(j, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
870 :
871 8610 : a[i - 1] +=
872 8610 : cp->e2[n] * term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
873 8610 : b[i - 1] +=
874 8610 : cp->n2[n] * term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
875 : }
876 : }
877 : }
878 :
879 80 : if (numactive <= m->n)
880 0 : return (MINTERR);
881 :
882 : /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
883 :
884 264 : for (int i = 2; i <= m->n; i++)
885 : {
886 520 : for (int j = 1; j < i; j++)
887 336 : M(i, j) = M(j, i);
888 : }
889 :
890 80 : return (solvemat(m, a, b, E, N));
891 : }
892 :
893 : /***************************************************************************/
894 : /*
895 : CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
896 :
897 : ORDER\TERM 1 2 3 4 5 6 7 8 9 10
898 : 1 e0n0 e1n0 e0n1
899 : 2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
900 : 3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
901 : */
902 : /***************************************************************************/
903 :
904 74844 : static double term(int nTerm, double e, double n)
905 : {
906 74844 : switch (nTerm)
907 : {
908 13464 : case 1:
909 13464 : return (1.0);
910 13464 : case 2:
911 13464 : return (e);
912 13464 : case 3:
913 13464 : return (n);
914 11484 : case 4:
915 11484 : return ((e * e));
916 11484 : case 5:
917 11484 : return ((e * n));
918 11484 : case 6:
919 11484 : return ((n * n));
920 0 : case 7:
921 0 : return ((e * e * e));
922 0 : case 8:
923 0 : return ((e * e * n));
924 0 : case 9:
925 0 : return ((e * n * n));
926 0 : case 10:
927 0 : return ((n * n * n));
928 : }
929 0 : return 0.0;
930 : }
931 :
932 : /***************************************************************************/
933 : /*
934 : SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
935 : GAUSSIAN ELIMINATION METHOD.
936 :
937 : | M11 M12 ... M1n | | E0 | | a0 |
938 : | M21 M22 ... M2n | | E1 | = | a1 |
939 : | . . . . | | . | | . |
940 : | Mn1 Mn2 ... Mnn | | En-1 | | an-1 |
941 :
942 : and
943 :
944 : | M11 M12 ... M1n | | N0 | | b0 |
945 : | M21 M22 ... M2n | | N1 | = | b1 |
946 : | . . . . | | . | | . |
947 : | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 |
948 : */
949 : /***************************************************************************/
950 :
951 104 : static int solvemat(struct MATRIX *m, double a[], double b[], double E[],
952 : double N[])
953 : {
954 440 : for (int i = 1; i <= m->n; i++)
955 : {
956 336 : int j = i;
957 :
958 : /* find row with largest magnitude value for pivot value */
959 :
960 336 : double pivot =
961 336 : M(i, j); /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
962 336 : int imark = i;
963 744 : for (int i2 = i + 1; i2 <= m->n; i2++)
964 : {
965 408 : if (fabs(M(i2, j)) > fabs(pivot))
966 : {
967 53 : pivot = M(i2, j);
968 53 : imark = i2;
969 : }
970 : }
971 :
972 : /* if the pivot is very small then the points are nearly co-linear */
973 : /* co-linear points result in an undefined matrix, and nearly */
974 : /* co-linear points results in a solution with rounding error */
975 :
976 336 : if (pivot == 0.0)
977 0 : return (MUNSOLVABLE);
978 :
979 : /* if row with highest pivot is not the current row, switch them */
980 :
981 336 : if (imark != i)
982 : {
983 245 : for (int j2 = 1; j2 <= m->n; j2++)
984 : {
985 198 : std::swap(M(imark, j2), M(i, j2));
986 : }
987 :
988 47 : std::swap(a[imark - 1], a[i - 1]);
989 47 : std::swap(b[imark - 1], b[i - 1]);
990 : }
991 :
992 : /* compute zeros above and below the pivot, and compute
993 : values for the rest of the row as well */
994 :
995 1488 : for (int i2 = 1; i2 <= m->n; i2++)
996 : {
997 1152 : if (i2 != i)
998 : {
999 816 : const double factor = M(i2, j) / pivot;
1000 2808 : for (int j2 = j; j2 <= m->n; j2++)
1001 1992 : M(i2, j2) -= factor * M(i, j2);
1002 816 : a[i2 - 1] -= factor * a[i - 1];
1003 816 : b[i2 - 1] -= factor * b[i - 1];
1004 : }
1005 : }
1006 : }
1007 :
1008 : /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
1009 : COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
1010 :
1011 440 : for (int i = 1; i <= m->n; i++)
1012 : {
1013 336 : E[i - 1] = a[i - 1] / M(i, i);
1014 336 : N[i - 1] = b[i - 1] / M(i, i);
1015 : }
1016 :
1017 104 : return (MSUCCESS);
1018 : }
1019 :
1020 : /***************************************************************************/
1021 : /*
1022 : DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
1023 : OUTLIER.
1024 :
1025 : THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
1026 : AND THE ALLOWED TOLERANCE:
1027 :
1028 : sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
1029 : lineAdj = b0 + b1*sample + b2*line + b3*line*sample
1030 :
1031 : WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
1032 :
1033 : [residualSample] = [A1][sampleCoefficients] - [b1]
1034 : [residualLine] = [A2][lineCoefficients] - [b2]
1035 :
1036 : sampleResidual^2 = sum( [residualSample]^2 )
1037 : lineResidual^2 = sum( [lineSample]^2 )
1038 :
1039 : residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
1040 :
1041 : THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
1042 : CONSIDERED THE WORST OUTLIER.
1043 :
1044 : IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
1045 : */
1046 : /***************************************************************************/
1047 3 : static int worst_outlier(struct Control_Points *cp, double x_mean,
1048 : double y_mean, int nOrder, double E[], double N[],
1049 : double dfTolerance)
1050 : {
1051 : // double dfSampleResidual = 0.0;
1052 : // double dfLineResidual = 0.0;
1053 : double *padfResiduals =
1054 3 : static_cast<double *>(CPLCalloc(sizeof(double), cp->count));
1055 :
1056 30 : for (int nI = 0; nI < cp->count; nI++)
1057 : {
1058 27 : double dfSampleRes = 0.0;
1059 27 : double dfLineRes = 0.0;
1060 27 : CRS_georef(cp->e1[nI] - x_mean, cp->n1[nI] - y_mean, &dfSampleRes,
1061 : &dfLineRes, E, N, nOrder);
1062 27 : dfSampleRes -= cp->e2[nI];
1063 27 : dfLineRes -= cp->n2[nI];
1064 : // dfSampleResidual += dfSampleRes*dfSampleRes;
1065 : // dfLineResidual += dfLineRes*dfLineRes;
1066 :
1067 27 : padfResiduals[nI] =
1068 27 : sqrt(dfSampleRes * dfSampleRes + dfLineRes * dfLineRes);
1069 : }
1070 :
1071 3 : int nIndex = -1;
1072 3 : double dfDifference = -1.0;
1073 30 : for (int nI = 0; nI < cp->count; nI++)
1074 : {
1075 27 : double dfCurrentDifference = padfResiduals[nI];
1076 27 : if (fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
1077 : {
1078 8 : dfCurrentDifference = 0.0;
1079 : }
1080 27 : if (dfCurrentDifference > dfDifference &&
1081 : dfCurrentDifference >= dfTolerance)
1082 : {
1083 5 : dfDifference = dfCurrentDifference;
1084 5 : nIndex = nI;
1085 : }
1086 : }
1087 3 : CPLFree(padfResiduals);
1088 3 : return nIndex;
1089 : }
1090 :
1091 : /***************************************************************************/
1092 : /*
1093 : REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
1094 : ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
1095 :
1096 : 1. WE CALCULATE THE COEFFICIENTS FOR ALL THE GCPS.
1097 : 2. THE GCP LIST WILL BE SCANNED TO DETERMINE THE WORST OUTLIER USING
1098 : THE CALCULATED COEFFICIENTS.
1099 : 3. THE WORST OUTLIER WILL BE REMOVED FROM THE GCP LIST.
1100 : 4. THE COEFFICIENTS WILL BE RECALCULATED WITHOUT THE WORST OUTLIER.
1101 : 5. STEP 1 TO 4 ARE EXECUTED UNTIL THE MINIMUM NUMBER OF GCPS WERE REACHED
1102 : OR IF NO GCP IS CONSIDERED AN OUTLIER WITH THE PASSED TOLERANCE.
1103 : */
1104 : /***************************************************************************/
1105 1 : static int remove_outliers(GCPTransformInfo *psInfo)
1106 : {
1107 1 : double *padfGeoX = nullptr;
1108 1 : double *padfGeoY = nullptr;
1109 1 : double *padfRasterX = nullptr;
1110 1 : double *padfRasterY = nullptr;
1111 1 : int *panStatus = nullptr;
1112 1 : int nCRSresult = 0;
1113 1 : int nGCPCount = 0;
1114 1 : int nMinimumGcps = 0;
1115 1 : int nReqOrder = 0;
1116 1 : double dfTolerance = 0;
1117 : struct Control_Points sPoints;
1118 :
1119 1 : double x1_sum = 0;
1120 1 : double y1_sum = 0;
1121 1 : double x2_sum = 0;
1122 1 : double y2_sum = 0;
1123 1 : memset(&sPoints, 0, sizeof(sPoints));
1124 :
1125 1 : nGCPCount = static_cast<int>(psInfo->asGCPs.size());
1126 1 : nMinimumGcps = psInfo->nMinimumGcps;
1127 1 : nReqOrder = psInfo->nOrder;
1128 1 : dfTolerance = psInfo->dfTolerance;
1129 :
1130 : try
1131 : {
1132 1 : padfGeoX = new double[nGCPCount];
1133 1 : padfGeoY = new double[nGCPCount];
1134 1 : padfRasterX = new double[nGCPCount];
1135 1 : padfRasterY = new double[nGCPCount];
1136 1 : panStatus = new int[nGCPCount];
1137 :
1138 11 : for (int nI = 0; nI < nGCPCount; nI++)
1139 : {
1140 10 : panStatus[nI] = 1;
1141 10 : padfGeoX[nI] = psInfo->asGCPs[nI].X();
1142 10 : padfGeoY[nI] = psInfo->asGCPs[nI].Y();
1143 10 : padfRasterX[nI] = psInfo->asGCPs[nI].Pixel();
1144 10 : padfRasterY[nI] = psInfo->asGCPs[nI].Line();
1145 10 : x1_sum += psInfo->asGCPs[nI].Pixel();
1146 10 : y1_sum += psInfo->asGCPs[nI].Line();
1147 10 : x2_sum += psInfo->asGCPs[nI].X();
1148 10 : y2_sum += psInfo->asGCPs[nI].Y();
1149 : }
1150 1 : psInfo->x1_mean = x1_sum / nGCPCount;
1151 1 : psInfo->y1_mean = y1_sum / nGCPCount;
1152 1 : psInfo->x2_mean = x2_sum / nGCPCount;
1153 1 : psInfo->y2_mean = y2_sum / nGCPCount;
1154 :
1155 1 : sPoints.count = nGCPCount;
1156 1 : sPoints.e1 = padfRasterX;
1157 1 : sPoints.n1 = padfRasterY;
1158 1 : sPoints.e2 = padfGeoX;
1159 1 : sPoints.n2 = padfGeoY;
1160 1 : sPoints.status = panStatus;
1161 :
1162 1 : nCRSresult = CRS_compute_georef_equations(
1163 1 : psInfo, &sPoints, psInfo->adfToGeoX, psInfo->adfToGeoY,
1164 1 : psInfo->adfFromGeoX, psInfo->adfFromGeoY, nReqOrder);
1165 :
1166 3 : while (sPoints.count > nMinimumGcps)
1167 : {
1168 6 : int nIndex = worst_outlier(
1169 : &sPoints, psInfo->x1_mean, psInfo->y1_mean, psInfo->nOrder,
1170 3 : psInfo->adfToGeoX, psInfo->adfToGeoY, dfTolerance);
1171 :
1172 : // If no outliers were detected, stop the GCP elimination
1173 3 : if (nIndex == -1)
1174 : {
1175 1 : break;
1176 : }
1177 :
1178 8 : for (int nI = nIndex; nI < sPoints.count - 1; nI++)
1179 : {
1180 6 : sPoints.e1[nI] = sPoints.e1[nI + 1];
1181 6 : sPoints.n1[nI] = sPoints.n1[nI + 1];
1182 6 : sPoints.e2[nI] = sPoints.e2[nI + 1];
1183 6 : sPoints.n2[nI] = sPoints.n2[nI + 1];
1184 6 : psInfo->asGCPs[nI].SetId(psInfo->asGCPs[nI + 1].Id());
1185 6 : psInfo->asGCPs[nI].SetInfo(psInfo->asGCPs[nI + 1].Info());
1186 : }
1187 :
1188 2 : sPoints.count = sPoints.count - 1;
1189 :
1190 2 : nCRSresult = CRS_compute_georef_equations(
1191 2 : psInfo, &sPoints, psInfo->adfToGeoX, psInfo->adfToGeoY,
1192 2 : psInfo->adfFromGeoX, psInfo->adfFromGeoY, nReqOrder);
1193 : }
1194 :
1195 9 : for (int nI = 0; nI < sPoints.count; nI++)
1196 : {
1197 8 : psInfo->asGCPs[nI].X() = sPoints.e2[nI];
1198 8 : psInfo->asGCPs[nI].Y() = sPoints.n2[nI];
1199 8 : psInfo->asGCPs[nI].Pixel() = sPoints.e1[nI];
1200 8 : psInfo->asGCPs[nI].Line() = sPoints.n1[nI];
1201 : }
1202 1 : psInfo->asGCPs.resize(sPoints.count);
1203 : }
1204 0 : catch (const std::exception &e)
1205 : {
1206 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
1207 0 : nCRSresult = MINTERR;
1208 : }
1209 1 : delete[] padfGeoX;
1210 1 : delete[] padfGeoY;
1211 1 : delete[] padfRasterX;
1212 1 : delete[] padfRasterY;
1213 1 : delete[] panStatus;
1214 :
1215 1 : return nCRSresult;
1216 : }
|
