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