Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Homography Transformer
4 : * Author: Nathan Olson
5 : *
6 : ******************************************************************************
7 : * Copyright (c) 2025, Nathan Olson <nathanmolson at gmail dot com>
8 : *
9 : * SPDX-License-Identifier: MIT
10 : ****************************************************************************/
11 :
12 : #include "cpl_port.h"
13 :
14 : #include <stdlib.h>
15 : #include <string.h>
16 : #include <algorithm>
17 :
18 : #include "cpl_atomic_ops.h"
19 : #include "cpl_error.h"
20 : #include "cpl_string.h"
21 : #include "gdal.h"
22 : #include "gdal_alg.h"
23 : #include "gdallinearsystem.h"
24 :
25 : CPL_C_START
26 : CPLXMLNode *GDALSerializeHomographyTransformer(void *pTransformArg);
27 : void *GDALDeserializeHomographyTransformer(CPLXMLNode *psTree);
28 : CPL_C_END
29 :
30 : struct HomographyTransformInfo
31 : {
32 : GDALTransformerInfo sTI{};
33 :
34 : double padfForward[9]{};
35 : double padfReverse[9]{};
36 :
37 : volatile int nRefCount{};
38 : };
39 :
40 : /************************************************************************/
41 : /* GDALCreateSimilarHomographyTransformer() */
42 : /************************************************************************/
43 :
44 0 : static void *GDALCreateSimilarHomographyTransformer(void *hTransformArg,
45 : double dfRatioX,
46 : double dfRatioY)
47 : {
48 0 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarHomographyTransformer",
49 : nullptr);
50 :
51 0 : HomographyTransformInfo *psInfo =
52 : static_cast<HomographyTransformInfo *>(hTransformArg);
53 :
54 0 : if (dfRatioX == 1.0 && dfRatioY == 1.0)
55 : {
56 : // We can just use a ref count, since using the source transformation
57 : // is thread-safe.
58 0 : CPLAtomicInc(&(psInfo->nRefCount));
59 : }
60 : else
61 : {
62 : double homography[9];
63 0 : for (int i = 0; i < 3; i++)
64 : {
65 0 : homography[3 * i + 1] = psInfo->padfForward[3 * i + 1] / dfRatioX;
66 0 : homography[3 * i + 2] = psInfo->padfForward[3 * i + 2] / dfRatioY;
67 0 : homography[3 * i] = psInfo->padfForward[3 * i];
68 : }
69 : psInfo = static_cast<HomographyTransformInfo *>(
70 0 : GDALCreateHomographyTransformer(homography));
71 : }
72 :
73 0 : return psInfo;
74 : }
75 :
76 : /************************************************************************/
77 : /* GDALCreateHomographyTransformer() */
78 : /************************************************************************/
79 :
80 : /**
81 : * Create Homography transformer from GCPs.
82 : *
83 : * Homography Transformers are serializable.
84 : *
85 : * @param adfHomography the forward homography.
86 : *
87 : * @return the transform argument or NULL if creation fails.
88 : */
89 :
90 29 : void *GDALCreateHomographyTransformer(double adfHomography[9])
91 : {
92 : /* -------------------------------------------------------------------- */
93 : /* Allocate transform info. */
94 : /* -------------------------------------------------------------------- */
95 29 : HomographyTransformInfo *psInfo = new HomographyTransformInfo();
96 :
97 29 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
98 : strlen(GDAL_GTI2_SIGNATURE));
99 29 : psInfo->sTI.pszClassName = "GDALHomographyTransformer";
100 29 : psInfo->sTI.pfnTransform = GDALHomographyTransform;
101 29 : psInfo->sTI.pfnCleanup = GDALDestroyHomographyTransformer;
102 29 : psInfo->sTI.pfnSerialize = GDALSerializeHomographyTransformer;
103 29 : psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarHomographyTransformer;
104 :
105 29 : psInfo->nRefCount = 1;
106 :
107 29 : memcpy(psInfo->padfForward, adfHomography, 9 * sizeof(double));
108 29 : if (GDALInvHomography(psInfo->padfForward, psInfo->padfReverse))
109 : {
110 29 : return psInfo;
111 : }
112 :
113 0 : GDALDestroyHomographyTransformer(psInfo);
114 0 : return nullptr;
115 : }
116 :
117 : /************************************************************************/
118 : /* GDALGCPsToHomography() */
119 : /************************************************************************/
120 :
121 : /**
122 : * \brief Generate Homography from GCPs.
123 : *
124 : * Given a set of GCPs perform least squares fit as a homography.
125 : *
126 : * A minimum of four GCPs are required to uniquely define a homography.
127 : * If there are less than four GCPs, GDALGCPsToGeoTransform() is used to
128 : * compute the transform.
129 : *
130 : * @param nGCPCount the number of GCPs being passed in.
131 : * @param pasGCPList the list of GCP structures.
132 : * @param padfHomography the nine double array in which the homography
133 : * will be returned.
134 : *
135 : * @return TRUE on success or FALSE if there aren't enough points to prepare a
136 : * homography, or pathological geometry is detected
137 : */
138 32 : int GDALGCPsToHomography(int nGCPCount, const GDAL_GCP *pasGCPList,
139 : double *padfHomography)
140 : {
141 32 : if (nGCPCount < 4)
142 : {
143 3 : padfHomography[6] = 1.0;
144 3 : padfHomography[7] = 0.0;
145 3 : padfHomography[8] = 0.0;
146 3 : return GDALGCPsToGeoTransform(nGCPCount, pasGCPList, padfHomography,
147 3 : FALSE);
148 : }
149 :
150 : /* -------------------------------------------------------------------- */
151 : /* Compute source and destination ranges so we can normalize */
152 : /* the values to make the least squares computation more stable. */
153 : /* -------------------------------------------------------------------- */
154 29 : double min_pixel = pasGCPList[0].dfGCPPixel;
155 29 : double max_pixel = pasGCPList[0].dfGCPPixel;
156 29 : double min_line = pasGCPList[0].dfGCPLine;
157 29 : double max_line = pasGCPList[0].dfGCPLine;
158 29 : double min_geox = pasGCPList[0].dfGCPX;
159 29 : double max_geox = pasGCPList[0].dfGCPX;
160 29 : double min_geoy = pasGCPList[0].dfGCPY;
161 29 : double max_geoy = pasGCPList[0].dfGCPY;
162 :
163 118 : for (int i = 1; i < nGCPCount; ++i)
164 : {
165 89 : min_pixel = std::min(min_pixel, pasGCPList[i].dfGCPPixel);
166 89 : max_pixel = std::max(max_pixel, pasGCPList[i].dfGCPPixel);
167 89 : min_line = std::min(min_line, pasGCPList[i].dfGCPLine);
168 89 : max_line = std::max(max_line, pasGCPList[i].dfGCPLine);
169 89 : min_geox = std::min(min_geox, pasGCPList[i].dfGCPX);
170 89 : max_geox = std::max(max_geox, pasGCPList[i].dfGCPX);
171 89 : min_geoy = std::min(min_geoy, pasGCPList[i].dfGCPY);
172 89 : max_geoy = std::max(max_geoy, pasGCPList[i].dfGCPY);
173 : }
174 :
175 29 : double EPS = 1.0e-12;
176 :
177 58 : if (std::abs(max_pixel - min_pixel) < EPS ||
178 58 : std::abs(max_line - min_line) < EPS ||
179 87 : std::abs(max_geox - min_geox) < EPS ||
180 29 : std::abs(max_geoy - min_geoy) < EPS)
181 : {
182 0 : return FALSE; // degenerate in at least one dimension.
183 : }
184 :
185 : double pl_normalize[9], geo_normalize[9];
186 :
187 29 : pl_normalize[0] = -min_pixel / (max_pixel - min_pixel);
188 29 : pl_normalize[1] = 1.0 / (max_pixel - min_pixel);
189 29 : pl_normalize[2] = 0.0;
190 29 : pl_normalize[3] = -min_line / (max_line - min_line);
191 29 : pl_normalize[4] = 0.0;
192 29 : pl_normalize[5] = 1.0 / (max_line - min_line);
193 29 : pl_normalize[6] = 1.0;
194 29 : pl_normalize[7] = 0.0;
195 29 : pl_normalize[8] = 0.0;
196 :
197 29 : geo_normalize[0] = -min_geox / (max_geox - min_geox);
198 29 : geo_normalize[1] = 1.0 / (max_geox - min_geox);
199 29 : geo_normalize[2] = 0.0;
200 29 : geo_normalize[3] = -min_geoy / (max_geoy - min_geoy);
201 29 : geo_normalize[4] = 0.0;
202 29 : geo_normalize[5] = 1.0 / (max_geoy - min_geoy);
203 29 : geo_normalize[6] = 1.0;
204 29 : geo_normalize[7] = 0.0;
205 29 : geo_normalize[8] = 0.0;
206 :
207 29 : double inv_geo_normalize[9] = {0.0};
208 29 : if (!GDALInvHomography(geo_normalize, inv_geo_normalize))
209 : {
210 0 : return FALSE;
211 : }
212 :
213 : /* -------------------------------------------------------------------- */
214 : /* Calculate the best fit homography following */
215 : /* https://www.cs.unc.edu/~ronisen/teaching/fall_2023/pdf_slides/ */
216 : /* lecture9_transformation.pdf */
217 : /* Since rank(AtA) = rank(8) = 8, append an additional equation */
218 : /* h_normalized[6] = 1 to fully define the solution. */
219 : /* -------------------------------------------------------------------- */
220 58 : GDALMatrix AtA(9, 9);
221 58 : GDALMatrix rhs(9, 1);
222 29 : rhs(6, 0) = 1;
223 29 : AtA(6, 6) = 1;
224 :
225 147 : for (int i = 0; i < nGCPCount; ++i)
226 : {
227 : double pixel, line, geox, geoy;
228 :
229 354 : if (!GDALApplyHomography(pl_normalize, pasGCPList[i].dfGCPPixel,
230 236 : pasGCPList[i].dfGCPLine, &pixel, &line) ||
231 118 : !GDALApplyHomography(geo_normalize, pasGCPList[i].dfGCPX,
232 118 : pasGCPList[i].dfGCPY, &geox, &geoy))
233 : {
234 0 : return FALSE;
235 : }
236 :
237 118 : double Ax[] = {1, pixel, line, 0, 0,
238 118 : 0, -geox, -geox * pixel, -geox * line};
239 118 : double Ay[] = {0, 0, 0, 1, pixel, line, -geoy, -geoy * pixel,
240 118 : -geoy * line};
241 : int j, k;
242 : // Populate the lower triangle of symmetric AtA matrix
243 1180 : for (j = 0; j < 9; j++)
244 : {
245 6372 : for (k = j; k < 9; k++)
246 : {
247 5310 : AtA(j, k) += Ax[j] * Ax[k] + Ay[j] * Ay[k];
248 : }
249 : }
250 : }
251 : // Populate the upper triangle of symmetric AtA matrix
252 290 : for (int j = 0; j < 9; j++)
253 : {
254 1305 : for (int k = 0; k < j; k++)
255 : {
256 1044 : AtA(j, k) = AtA(k, j);
257 : }
258 : }
259 :
260 58 : GDALMatrix h_normalized(9, 1);
261 29 : if (!GDALLinearSystemSolve(AtA, rhs, h_normalized))
262 : {
263 0 : return FALSE;
264 : }
265 29 : if (abs(h_normalized(6, 0)) < 1.0e-15)
266 : {
267 0 : return FALSE;
268 : }
269 :
270 : /* -------------------------------------------------------------------- */
271 : /* Check that the homography maps the unit square to a convex */
272 : /* quadrilateral. */
273 : /* -------------------------------------------------------------------- */
274 : // First, use the normalized homography to make the corners of the unit
275 : // square to normalized geo coordinates
276 29 : double x[4] = {0, 1, 1, 0};
277 29 : double y[4] = {0, 0, 1, 1};
278 145 : for (int i = 0; i < 4; i++)
279 : {
280 116 : if (!GDALApplyHomography(h_normalized.data(), x[i], y[i], &x[i], &y[i]))
281 : {
282 0 : return FALSE;
283 : }
284 : }
285 : // Next, compute the vector from the top-left corner to each corner
286 145 : for (int i = 3; i >= 0; i--)
287 : {
288 116 : x[i] -= x[0];
289 116 : y[i] -= y[0];
290 : }
291 : // Finally, check that "v2" (<x[2], y[2]>, the vector from top-left to
292 : // bottom-right corner) is between v1 and v3, by checking that the
293 : // vector cross product (v1 x v2) has the same sign as (v2 x v3)
294 29 : double cross12 = x[1] * y[2] - x[2] * y[1];
295 29 : double cross23 = x[2] * y[3] - x[3] * y[2];
296 29 : if (cross12 * cross23 <= 0.0)
297 : {
298 2 : return FALSE;
299 : }
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Compose the resulting transformation with the normalization */
303 : /* homographies. */
304 : /* -------------------------------------------------------------------- */
305 27 : double h1p2[9] = {0.0};
306 :
307 27 : GDALComposeHomographies(pl_normalize, h_normalized.data(), h1p2);
308 27 : GDALComposeHomographies(h1p2, inv_geo_normalize, padfHomography);
309 :
310 27 : return TRUE;
311 : }
312 :
313 : /************************************************************************/
314 : /* GDALComposeHomographies() */
315 : /************************************************************************/
316 :
317 : /**
318 : * \brief Compose two homographies.
319 : *
320 : * The resulting homography is the equivalent to padfH1 and then padfH2
321 : * being applied to a point.
322 : *
323 : * @param padfH1 the first homography, nine values.
324 : * @param padfH2 the second homography, nine values.
325 : * @param padfHOut the output homography, nine values, may safely be the same
326 : * array as padfH1 or padfH2.
327 : */
328 :
329 54 : void GDALComposeHomographies(const double *padfH1, const double *padfH2,
330 : double *padfHOut)
331 :
332 : {
333 54 : double hwrk[9] = {0.0};
334 :
335 54 : hwrk[1] =
336 54 : padfH2[1] * padfH1[1] + padfH2[2] * padfH1[4] + padfH2[0] * padfH1[7];
337 54 : hwrk[2] =
338 54 : padfH2[1] * padfH1[2] + padfH2[2] * padfH1[5] + padfH2[0] * padfH1[8];
339 54 : hwrk[0] =
340 54 : padfH2[1] * padfH1[0] + padfH2[2] * padfH1[3] + padfH2[0] * padfH1[6];
341 :
342 54 : hwrk[4] =
343 54 : padfH2[4] * padfH1[1] + padfH2[5] * padfH1[4] + padfH2[3] * padfH1[7];
344 54 : hwrk[5] =
345 54 : padfH2[4] * padfH1[2] + padfH2[5] * padfH1[5] + padfH2[3] * padfH1[8];
346 54 : hwrk[3] =
347 54 : padfH2[4] * padfH1[0] + padfH2[5] * padfH1[3] + padfH2[3] * padfH1[6];
348 :
349 54 : hwrk[7] =
350 54 : padfH2[7] * padfH1[1] + padfH2[8] * padfH1[4] + padfH2[6] * padfH1[7];
351 54 : hwrk[8] =
352 54 : padfH2[7] * padfH1[2] + padfH2[8] * padfH1[5] + padfH2[6] * padfH1[8];
353 54 : hwrk[6] =
354 54 : padfH2[7] * padfH1[0] + padfH2[8] * padfH1[3] + padfH2[6] * padfH1[6];
355 54 : memcpy(padfHOut, hwrk, sizeof(hwrk));
356 54 : }
357 :
358 : /************************************************************************/
359 : /* GDALApplyHomography() */
360 : /************************************************************************/
361 :
362 : /**
363 : * Apply Homography to x/y coordinate.
364 : *
365 : * Applies the following computation, converting a (pixel, line) coordinate
366 : * into a georeferenced (geo_x, geo_y) location.
367 : * \code{.c}
368 : * *pdfGeoX = (padfHomography[0] + dfPixel * padfHomography[1]
369 : * + dfLine * padfHomography[2]) /
370 : * (padfHomography[6] + dfPixel * padfHomography[7]
371 : * + dfLine * padfHomography[8]);
372 : * *pdfGeoY = (padfHomography[3] + dfPixel * padfHomography[4]
373 : * + dfLine * padfHomography[5]) /
374 : * (padfHomography[6] + dfPixel * padfHomography[7]
375 : * + dfLine * padfHomography[8]);
376 : * \endcode
377 : *
378 : * @param padfHomography Nine coefficient Homography to apply.
379 : * @param dfPixel Input pixel position.
380 : * @param dfLine Input line position.
381 : * @param pdfGeoX output location where geo_x (easting/longitude)
382 : * location is placed.
383 : * @param pdfGeoY output location where geo_y (northing/latitude)
384 : * location is placed.
385 : *
386 : * @return TRUE on success or FALSE if failure.
387 : */
388 :
389 362 : int GDALApplyHomography(const double *padfHomography, double dfPixel,
390 : double dfLine, double *pdfGeoX, double *pdfGeoY)
391 : {
392 362 : double w = padfHomography[6] + dfPixel * padfHomography[7] +
393 362 : dfLine * padfHomography[8];
394 362 : if (abs(w) < 1.0e-15)
395 : {
396 0 : return FALSE;
397 : }
398 362 : double wx = padfHomography[0] + dfPixel * padfHomography[1] +
399 362 : dfLine * padfHomography[2];
400 362 : double wy = padfHomography[3] + dfPixel * padfHomography[4] +
401 362 : dfLine * padfHomography[5];
402 362 : *pdfGeoX = wx / w;
403 362 : *pdfGeoY = wy / w;
404 362 : return TRUE;
405 : }
406 :
407 : /************************************************************************/
408 : /* GDALInvHomography() */
409 : /************************************************************************/
410 :
411 : /**
412 : * Invert Homography.
413 : *
414 : * This function will invert a standard 3x3 set of Homography coefficients.
415 : * This converts the equation from being pixel to geo to being geo to pixel.
416 : *
417 : * @param padfHIn Input homography (nine doubles - unaltered).
418 : * @param padfHOut Output homography (nine doubles - updated).
419 : *
420 : * @return TRUE on success or FALSE if the equation is uninvertable.
421 : */
422 :
423 65 : int GDALInvHomography(const double *padfHIn, double *padfHOut)
424 :
425 : {
426 : // Special case - no rotation - to avoid computing determinate
427 : // and potential precision issues.
428 65 : if (padfHIn[2] == 0.0 && padfHIn[4] == 0.0 && padfHIn[1] != 0.0 &&
429 59 : padfHIn[5] != 0.0 && padfHIn[7] == 0.0 && padfHIn[8] == 0.0 &&
430 59 : padfHIn[6] != 0.0)
431 : {
432 59 : padfHOut[0] = -padfHIn[0] / padfHIn[1] / padfHIn[6];
433 59 : padfHOut[1] = 1.0 / padfHIn[1];
434 59 : padfHOut[2] = 0.0;
435 59 : padfHOut[3] = -padfHIn[3] / padfHIn[5] / padfHIn[6];
436 59 : padfHOut[4] = 0.0;
437 59 : padfHOut[5] = 1.0 / padfHIn[5];
438 59 : padfHOut[6] = 1.0 / padfHIn[6];
439 59 : padfHOut[7] = 0.0;
440 59 : padfHOut[8] = 0.0;
441 59 : return 1;
442 : }
443 :
444 : // Compute determinate.
445 :
446 6 : const double det = padfHIn[1] * padfHIn[5] * padfHIn[6] -
447 6 : padfHIn[2] * padfHIn[4] * padfHIn[6] +
448 6 : padfHIn[2] * padfHIn[3] * padfHIn[7] -
449 6 : padfHIn[0] * padfHIn[5] * padfHIn[7] +
450 6 : padfHIn[0] * padfHIn[4] * padfHIn[8] -
451 6 : padfHIn[1] * padfHIn[3] * padfHIn[8];
452 : const double magnitude =
453 12 : std::max(std::max(fabs(padfHIn[1]), fabs(padfHIn[2])),
454 6 : std::max(fabs(padfHIn[4]), fabs(padfHIn[5])));
455 :
456 6 : if (fabs(det) <= 1e-10 * magnitude * magnitude)
457 3 : return 0;
458 :
459 3 : const double inv_det = 1.0 / det;
460 :
461 : // Compute adjoint, and divide by determinant.
462 :
463 3 : padfHOut[1] = (padfHIn[5] * padfHIn[6] - padfHIn[3] * padfHIn[8]) * inv_det;
464 3 : padfHOut[4] = (padfHIn[3] * padfHIn[7] - padfHIn[4] * padfHIn[6]) * inv_det;
465 3 : padfHOut[7] = (padfHIn[4] * padfHIn[8] - padfHIn[5] * padfHIn[7]) * inv_det;
466 :
467 3 : padfHOut[2] = (padfHIn[0] * padfHIn[8] - padfHIn[2] * padfHIn[6]) * inv_det;
468 3 : padfHOut[5] = (padfHIn[1] * padfHIn[6] - padfHIn[0] * padfHIn[7]) * inv_det;
469 3 : padfHOut[8] = (padfHIn[2] * padfHIn[7] - padfHIn[1] * padfHIn[8]) * inv_det;
470 :
471 3 : padfHOut[0] = (padfHIn[2] * padfHIn[3] - padfHIn[0] * padfHIn[5]) * inv_det;
472 3 : padfHOut[3] = (padfHIn[0] * padfHIn[4] - padfHIn[1] * padfHIn[3]) * inv_det;
473 3 : padfHOut[6] = (padfHIn[1] * padfHIn[5] - padfHIn[2] * padfHIn[4]) * inv_det;
474 :
475 3 : return 1;
476 : }
477 :
478 : /************************************************************************/
479 : /* GDALCreateHomographyTransformerFromGCPs() */
480 : /************************************************************************/
481 :
482 : /**
483 : * Create Homography transformer from GCPs.
484 : *
485 : * Homography Transformers are serializable.
486 : *
487 : * @param nGCPCount the number of GCPs in pasGCPList.
488 : * @param pasGCPList an array of GCPs to be used as input.
489 : *
490 : * @return the transform argument or NULL if creation fails.
491 : */
492 :
493 22 : void *GDALCreateHomographyTransformerFromGCPs(int nGCPCount,
494 : const GDAL_GCP *pasGCPList)
495 : {
496 : double adfHomography[9];
497 :
498 22 : if (GDALGCPsToHomography(nGCPCount, pasGCPList, adfHomography))
499 : {
500 22 : return GDALCreateHomographyTransformer(adfHomography);
501 : }
502 0 : return nullptr;
503 : }
504 :
505 : /************************************************************************/
506 : /* GDALDestroyHomographyTransformer() */
507 : /************************************************************************/
508 :
509 : /**
510 : * Destroy Homography transformer.
511 : *
512 : * This function is used to destroy information about a homography
513 : * transformation created with GDALCreateHomographyTransformer().
514 : *
515 : * @param pTransformArg the transform arg previously returned by
516 : * GDALCreateHomographyTransformer().
517 : */
518 :
519 29 : void GDALDestroyHomographyTransformer(void *pTransformArg)
520 :
521 : {
522 29 : if (pTransformArg == nullptr)
523 0 : return;
524 :
525 29 : HomographyTransformInfo *psInfo =
526 : static_cast<HomographyTransformInfo *>(pTransformArg);
527 :
528 29 : if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
529 : {
530 29 : delete psInfo;
531 : }
532 : }
533 :
534 : /************************************************************************/
535 : /* GDALHomographyTransform() */
536 : /************************************************************************/
537 :
538 : /**
539 : * Transforms point based on homography.
540 : *
541 : * This function matches the GDALTransformerFunc signature, and can be
542 : * used to transform one or more points from pixel/line coordinates to
543 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
544 : *
545 : * @param pTransformArg return value from GDALCreateHomographyTransformer().
546 : * @param bDstToSrc TRUE if transformation is from the destination
547 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
548 : * from pixel/line to georeferenced coordinates.
549 : * @param nPointCount the number of values in the x, y and z arrays.
550 : * @param x array containing the X values to be transformed.
551 : * @param y array containing the Y values to be transformed.
552 : * @param z array containing the Z values to be transformed.
553 : * @param panSuccess array in which a flag indicating success (TRUE) or
554 : * failure (FALSE) of the transformation are placed.
555 : *
556 : * @return TRUE if all points have been successfully transformed.
557 : */
558 :
559 17543 : int GDALHomographyTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
560 : double *x, double *y, CPL_UNUSED double *z,
561 : int *panSuccess)
562 : {
563 17543 : VALIDATE_POINTER1(pTransformArg, "GDALHomographyTransform", 0);
564 :
565 17543 : HomographyTransformInfo *psInfo =
566 : static_cast<HomographyTransformInfo *>(pTransformArg);
567 :
568 17543 : double *homography = bDstToSrc ? psInfo->padfReverse : psInfo->padfForward;
569 17543 : int ret = TRUE;
570 49853 : for (int i = 0; i < nPointCount; i++)
571 : {
572 32310 : double w = homography[6] + x[i] * homography[7] + y[i] * homography[8];
573 32310 : if (abs(w) < 1.0e-15)
574 : {
575 0 : panSuccess[i] = FALSE;
576 0 : ret = FALSE;
577 : }
578 : else
579 : {
580 32310 : double wx =
581 32310 : homography[0] + x[i] * homography[1] + y[i] * homography[2];
582 32310 : double wy =
583 32310 : homography[3] + x[i] * homography[4] + y[i] * homography[5];
584 32310 : x[i] = wx / w;
585 32310 : y[i] = wy / w;
586 32310 : panSuccess[i] = TRUE;
587 : }
588 : }
589 :
590 17543 : return ret;
591 : }
592 :
593 : /************************************************************************/
594 : /* GDALSerializeHomographyTransformer() */
595 : /************************************************************************/
596 :
597 8 : CPLXMLNode *GDALSerializeHomographyTransformer(void *pTransformArg)
598 :
599 : {
600 8 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeHomographyTransformer",
601 : nullptr);
602 :
603 8 : HomographyTransformInfo *psInfo =
604 : static_cast<HomographyTransformInfo *>(pTransformArg);
605 :
606 : CPLXMLNode *psTree =
607 8 : CPLCreateXMLNode(nullptr, CXT_Element, "HomographyTransformer");
608 :
609 : /* -------------------------------------------------------------------- */
610 : /* Attach Homography. */
611 : /* -------------------------------------------------------------------- */
612 8 : char szWork[300] = {};
613 :
614 8 : CPLsnprintf(
615 : szWork, sizeof(szWork),
616 : "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
617 : psInfo->padfForward[0], psInfo->padfForward[1], psInfo->padfForward[2],
618 : psInfo->padfForward[3], psInfo->padfForward[4], psInfo->padfForward[5],
619 : psInfo->padfForward[6], psInfo->padfForward[7], psInfo->padfForward[8]);
620 8 : CPLCreateXMLElementAndValue(psTree, "Homography", szWork);
621 :
622 8 : return psTree;
623 : }
624 :
625 : /************************************************************************/
626 : /* GDALDeserializeHomography() */
627 : /************************************************************************/
628 :
629 7 : static void GDALDeserializeHomography(const char *pszH, double adfHomography[9])
630 : {
631 7 : CPLsscanf(pszH, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf", adfHomography + 0,
632 : adfHomography + 1, adfHomography + 2, adfHomography + 3,
633 : adfHomography + 4, adfHomography + 5, adfHomography + 6,
634 : adfHomography + 7, adfHomography + 8);
635 7 : }
636 :
637 : /************************************************************************/
638 : /* GDALDeserializeHomographyTransformer() */
639 : /************************************************************************/
640 :
641 7 : void *GDALDeserializeHomographyTransformer(CPLXMLNode *psTree)
642 :
643 : {
644 : /* -------------------------------------------------------------------- */
645 : /* Homography */
646 : /* -------------------------------------------------------------------- */
647 : double padfForward[9];
648 7 : if (CPLGetXMLNode(psTree, "Homography") != nullptr)
649 : {
650 7 : GDALDeserializeHomography(CPLGetXMLValue(psTree, "Homography", ""),
651 : padfForward);
652 :
653 : /* -------------------------------------------------------------------- */
654 : /* Generate transformation. */
655 : /* -------------------------------------------------------------------- */
656 7 : void *pResult = GDALCreateHomographyTransformer(padfForward);
657 :
658 7 : return pResult;
659 : }
660 0 : return nullptr;
661 : }
|
