• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OSGeo / gdal / 13836648005

13 Mar 2025 02:09PM UTC coverage: 70.436% (-0.01%) from 70.446%
13836648005

push

github

web-flow
New Transform type: Homography (#11949)

Add new transform type, Homography.
Add functions to compute homography from a list of GCPs.
Add functions to serialize and deserialize a homography
Automatically select homography transfrom when there are 4 or 5 GCPs present.

Fixes #11940

231 of 274 new or added lines in 2 files covered. (84.31%)

16257 existing lines in 42 files now uncovered.

553736 of 786159 relevant lines covered (70.44%)

221595.72 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

89.18
/alg/gdal_homography.cpp
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

NEW
44
static void *GDALCreateSimilarHomographyTransformer(void *hTransformArg,
×
45
                                                    double dfRatioX,
46
                                                    double dfRatioY)
47
{
NEW
48
    VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarHomographyTransformer",
×
49
                      nullptr);
50

NEW
51
    HomographyTransformInfo *psInfo =
×
52
        static_cast<HomographyTransformInfo *>(hTransformArg);
53

NEW
54
    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.
NEW
58
        CPLAtomicInc(&(psInfo->nRefCount));
×
59
    }
60
    else
61
    {
62
        double homography[9];
NEW
63
        for (int i = 0; i < 3; i++)
×
64
        {
NEW
65
            homography[3 * i + 1] = psInfo->padfForward[3 * i + 1] / dfRatioX;
×
NEW
66
            homography[3 * i + 2] = psInfo->padfForward[3 * i + 2] / dfRatioY;
×
NEW
67
            homography[3 * i] = psInfo->padfForward[3 * i];
×
68
        }
69
        psInfo = static_cast<HomographyTransformInfo *>(
NEW
70
            GDALCreateHomographyTransformer(homography));
×
71
    }
72

NEW
73
    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
void *GDALCreateHomographyTransformer(double adfHomography[9])
29✔
91
{
92
    /* -------------------------------------------------------------------- */
93
    /*      Allocate transform info.                                        */
94
    /* -------------------------------------------------------------------- */
95
    HomographyTransformInfo *psInfo = new HomographyTransformInfo();
29✔
96

97
    memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
29✔
98
           strlen(GDAL_GTI2_SIGNATURE));
99
    psInfo->sTI.pszClassName = "GDALHomographyTransformer";
29✔
100
    psInfo->sTI.pfnTransform = GDALHomographyTransform;
29✔
101
    psInfo->sTI.pfnCleanup = GDALDestroyHomographyTransformer;
29✔
102
    psInfo->sTI.pfnSerialize = GDALSerializeHomographyTransformer;
29✔
103
    psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarHomographyTransformer;
29✔
104

105
    psInfo->nRefCount = 1;
29✔
106

107
    memcpy(psInfo->padfForward, adfHomography, 9 * sizeof(double));
29✔
108
    if (GDALInvHomography(psInfo->padfForward, psInfo->padfReverse))
29✔
109
    {
110
        return psInfo;
29✔
111
    }
112

NEW
113
    GDALDestroyHomographyTransformer(psInfo);
×
NEW
114
    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
int GDALGCPsToHomography(int nGCPCount, const GDAL_GCP *pasGCPList,
32✔
139
                         double *padfHomography)
140
{
141
    if (nGCPCount < 4)
32✔
142
    {
143
        padfHomography[6] = 1.0;
3✔
144
        padfHomography[7] = 0.0;
3✔
145
        padfHomography[8] = 0.0;
3✔
146
        return GDALGCPsToGeoTransform(nGCPCount, pasGCPList, padfHomography,
3✔
147
                                      FALSE);
3✔
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
    double min_pixel = pasGCPList[0].dfGCPPixel;
29✔
155
    double max_pixel = pasGCPList[0].dfGCPPixel;
29✔
156
    double min_line = pasGCPList[0].dfGCPLine;
29✔
157
    double max_line = pasGCPList[0].dfGCPLine;
29✔
158
    double min_geox = pasGCPList[0].dfGCPX;
29✔
159
    double max_geox = pasGCPList[0].dfGCPX;
29✔
160
    double min_geoy = pasGCPList[0].dfGCPY;
29✔
161
    double max_geoy = pasGCPList[0].dfGCPY;
29✔
162

163
    for (int i = 1; i < nGCPCount; ++i)
118✔
164
    {
165
        min_pixel = std::min(min_pixel, pasGCPList[i].dfGCPPixel);
89✔
166
        max_pixel = std::max(max_pixel, pasGCPList[i].dfGCPPixel);
89✔
167
        min_line = std::min(min_line, pasGCPList[i].dfGCPLine);
89✔
168
        max_line = std::max(max_line, pasGCPList[i].dfGCPLine);
89✔
169
        min_geox = std::min(min_geox, pasGCPList[i].dfGCPX);
89✔
170
        max_geox = std::max(max_geox, pasGCPList[i].dfGCPX);
89✔
171
        min_geoy = std::min(min_geoy, pasGCPList[i].dfGCPY);
89✔
172
        max_geoy = std::max(max_geoy, pasGCPList[i].dfGCPY);
89✔
173
    }
174

175
    double EPS = 1.0e-12;
29✔
176

177
    if (std::abs(max_pixel - min_pixel) < EPS ||
58✔
178
        std::abs(max_line - min_line) < EPS ||
58✔
179
        std::abs(max_geox - min_geox) < EPS ||
87✔
180
        std::abs(max_geoy - min_geoy) < EPS)
29✔
181
    {
NEW
182
        return FALSE;  // degenerate in at least one dimension.
×
183
    }
184

185
    double pl_normalize[9], geo_normalize[9];
186

187
    pl_normalize[0] = -min_pixel / (max_pixel - min_pixel);
29✔
188
    pl_normalize[1] = 1.0 / (max_pixel - min_pixel);
29✔
189
    pl_normalize[2] = 0.0;
29✔
190
    pl_normalize[3] = -min_line / (max_line - min_line);
29✔
191
    pl_normalize[4] = 0.0;
29✔
192
    pl_normalize[5] = 1.0 / (max_line - min_line);
29✔
193
    pl_normalize[6] = 1.0;
29✔
194
    pl_normalize[7] = 0.0;
29✔
195
    pl_normalize[8] = 0.0;
29✔
196

197
    geo_normalize[0] = -min_geox / (max_geox - min_geox);
29✔
198
    geo_normalize[1] = 1.0 / (max_geox - min_geox);
29✔
199
    geo_normalize[2] = 0.0;
29✔
200
    geo_normalize[3] = -min_geoy / (max_geoy - min_geoy);
29✔
201
    geo_normalize[4] = 0.0;
29✔
202
    geo_normalize[5] = 1.0 / (max_geoy - min_geoy);
29✔
203
    geo_normalize[6] = 1.0;
29✔
204
    geo_normalize[7] = 0.0;
29✔
205
    geo_normalize[8] = 0.0;
29✔
206

207
    double inv_geo_normalize[9] = {0.0};
29✔
208
    if (!GDALInvHomography(geo_normalize, inv_geo_normalize))
29✔
209
    {
NEW
210
        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
    GDALMatrix AtA(9, 9);
58✔
221
    GDALMatrix rhs(9, 1);
58✔
222
    rhs(6, 0) = 1;
29✔
223
    AtA(6, 6) = 1;
29✔
224

225
    for (int i = 0; i < nGCPCount; ++i)
147✔
226
    {
227
        double pixel, line, geox, geoy;
228

229
        if (!GDALApplyHomography(pl_normalize, pasGCPList[i].dfGCPPixel,
354✔
230
                                 pasGCPList[i].dfGCPLine, &pixel, &line) ||
236✔
231
            !GDALApplyHomography(geo_normalize, pasGCPList[i].dfGCPX,
118✔
232
                                 pasGCPList[i].dfGCPY, &geox, &geoy))
118✔
233
        {
NEW
234
            return FALSE;
×
235
        }
236

237
        double Ax[] = {1, pixel, line,          0,           0,
118✔
238
                       0, -geox, -geox * pixel, -geox * line};
118✔
239
        double Ay[] = {0,           0, 0, 1, pixel, line, -geoy, -geoy * pixel,
118✔
240
                       -geoy * line};
118✔
241
        int j, k;
242
        // Populate the lower triangle of symmetric AtA matrix
243
        for (j = 0; j < 9; j++)
1,180✔
244
        {
245
            for (k = j; k < 9; k++)
6,372✔
246
            {
247
                AtA(j, k) += Ax[j] * Ax[k] + Ay[j] * Ay[k];
5,310✔
248
            }
249
        }
250
    }
251
    // Populate the upper triangle of symmetric AtA matrix
252
    for (int j = 0; j < 9; j++)
290✔
253
    {
254
        for (int k = 0; k < j; k++)
1,305✔
255
        {
256
            AtA(j, k) = AtA(k, j);
1,044✔
257
        }
258
    }
259

260
    GDALMatrix h_normalized(9, 1);
58✔
261
    if (!GDALLinearSystemSolve(AtA, rhs, h_normalized))
29✔
262
    {
NEW
263
        return FALSE;
×
264
    }
265
    if (abs(h_normalized(6, 0)) < 1.0e-15)
29✔
266
    {
NEW
267
        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
    double x[4] = {0, 1, 1, 0};
29✔
277
    double y[4] = {0, 0, 1, 1};
29✔
278
    for (int i = 0; i < 4; i++)
145✔
279
    {
280
        if (!GDALApplyHomography(h_normalized.data(), x[i], y[i], &x[i], &y[i]))
116✔
281
        {
NEW
282
            return FALSE;
×
283
        }
284
    }
285
    // Next, compute the vector from the top-left corner to each corner
286
    for (int i = 3; i >= 0; i--)
145✔
287
    {
288
        x[i] -= x[0];
116✔
289
        y[i] -= y[0];
116✔
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
    double cross12 = x[1] * y[2] - x[2] * y[1];
29✔
295
    double cross23 = x[2] * y[3] - x[3] * y[2];
29✔
296
    if (cross12 * cross23 <= 0.0)
29✔
297
    {
298
        return FALSE;
2✔
299
    }
300

301
    /* -------------------------------------------------------------------- */
302
    /*      Compose the resulting transformation with the normalization     */
303
    /*      homographies.                                                   */
304
    /* -------------------------------------------------------------------- */
305
    double h1p2[9] = {0.0};
27✔
306

307
    GDALComposeHomographies(pl_normalize, h_normalized.data(), h1p2);
27✔
308
    GDALComposeHomographies(h1p2, inv_geo_normalize, padfHomography);
27✔
309

310
    return TRUE;
27✔
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
void GDALComposeHomographies(const double *padfH1, const double *padfH2,
54✔
330
                             double *padfHOut)
331

332
{
333
    double hwrk[9] = {0.0};
54✔
334

335
    hwrk[1] =
54✔
336
        padfH2[1] * padfH1[1] + padfH2[2] * padfH1[4] + padfH2[0] * padfH1[7];
54✔
337
    hwrk[2] =
54✔
338
        padfH2[1] * padfH1[2] + padfH2[2] * padfH1[5] + padfH2[0] * padfH1[8];
54✔
339
    hwrk[0] =
54✔
340
        padfH2[1] * padfH1[0] + padfH2[2] * padfH1[3] + padfH2[0] * padfH1[6];
54✔
341

342
    hwrk[4] =
54✔
343
        padfH2[4] * padfH1[1] + padfH2[5] * padfH1[4] + padfH2[3] * padfH1[7];
54✔
344
    hwrk[5] =
54✔
345
        padfH2[4] * padfH1[2] + padfH2[5] * padfH1[5] + padfH2[3] * padfH1[8];
54✔
346
    hwrk[3] =
54✔
347
        padfH2[4] * padfH1[0] + padfH2[5] * padfH1[3] + padfH2[3] * padfH1[6];
54✔
348

349
    hwrk[7] =
54✔
350
        padfH2[7] * padfH1[1] + padfH2[8] * padfH1[4] + padfH2[6] * padfH1[7];
54✔
351
    hwrk[8] =
54✔
352
        padfH2[7] * padfH1[2] + padfH2[8] * padfH1[5] + padfH2[6] * padfH1[8];
54✔
353
    hwrk[6] =
54✔
354
        padfH2[7] * padfH1[0] + padfH2[8] * padfH1[3] + padfH2[6] * padfH1[6];
54✔
355
    memcpy(padfHOut, hwrk, sizeof(hwrk));
54✔
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
int GDALApplyHomography(const double *padfHomography, double dfPixel,
362✔
390
                        double dfLine, double *pdfGeoX, double *pdfGeoY)
391
{
392
    double w = padfHomography[6] + dfPixel * padfHomography[7] +
362✔
393
               dfLine * padfHomography[8];
362✔
394
    if (abs(w) < 1.0e-15)
362✔
395
    {
NEW
396
        return FALSE;
×
397
    }
398
    double wx = padfHomography[0] + dfPixel * padfHomography[1] +
362✔
399
                dfLine * padfHomography[2];
362✔
400
    double wy = padfHomography[3] + dfPixel * padfHomography[4] +
362✔
401
                dfLine * padfHomography[5];
362✔
402
    *pdfGeoX = wx / w;
362✔
403
    *pdfGeoY = wy / w;
362✔
404
    return TRUE;
362✔
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
int GDALInvHomography(const double *padfHIn, double *padfHOut)
65✔
424

425
{
426
    // Special case - no rotation - to avoid computing determinate
427
    // and potential precision issues.
428
    if (padfHIn[2] == 0.0 && padfHIn[4] == 0.0 && padfHIn[1] != 0.0 &&
65✔
429
        padfHIn[5] != 0.0 && padfHIn[7] == 0.0 && padfHIn[8] == 0.0 &&
59✔
430
        padfHIn[6] != 0.0)
59✔
431
    {
432
        padfHOut[0] = -padfHIn[0] / padfHIn[1] / padfHIn[6];
59✔
433
        padfHOut[1] = 1.0 / padfHIn[1];
59✔
434
        padfHOut[2] = 0.0;
59✔
435
        padfHOut[3] = -padfHIn[3] / padfHIn[5] / padfHIn[6];
59✔
436
        padfHOut[4] = 0.0;
59✔
437
        padfHOut[5] = 1.0 / padfHIn[5];
59✔
438
        padfHOut[6] = 1.0 / padfHIn[6];
59✔
439
        padfHOut[7] = 0.0;
59✔
440
        padfHOut[8] = 0.0;
59✔
441
        return 1;
59✔
442
    }
443

444
    // Compute determinate.
445

446
    const double det = padfHIn[1] * padfHIn[5] * padfHIn[6] -
6✔
447
                       padfHIn[2] * padfHIn[4] * padfHIn[6] +
6✔
448
                       padfHIn[2] * padfHIn[3] * padfHIn[7] -
6✔
449
                       padfHIn[0] * padfHIn[5] * padfHIn[7] +
6✔
450
                       padfHIn[0] * padfHIn[4] * padfHIn[8] -
6✔
451
                       padfHIn[1] * padfHIn[3] * padfHIn[8];
6✔
452
    const double magnitude =
453
        std::max(std::max(fabs(padfHIn[1]), fabs(padfHIn[2])),
12✔
454
                 std::max(fabs(padfHIn[4]), fabs(padfHIn[5])));
6✔
455

456
    if (fabs(det) <= 1e-10 * magnitude * magnitude)
6✔
457
        return 0;
3✔
458

459
    const double inv_det = 1.0 / det;
3✔
460

461
    // Compute adjoint, and divide by determinant.
462

463
    padfHOut[1] = (padfHIn[5] * padfHIn[6] - padfHIn[3] * padfHIn[8]) * inv_det;
3✔
464
    padfHOut[4] = (padfHIn[3] * padfHIn[7] - padfHIn[4] * padfHIn[6]) * inv_det;
3✔
465
    padfHOut[7] = (padfHIn[4] * padfHIn[8] - padfHIn[5] * padfHIn[7]) * inv_det;
3✔
466

467
    padfHOut[2] = (padfHIn[0] * padfHIn[8] - padfHIn[2] * padfHIn[6]) * inv_det;
3✔
468
    padfHOut[5] = (padfHIn[1] * padfHIn[6] - padfHIn[0] * padfHIn[7]) * inv_det;
3✔
469
    padfHOut[8] = (padfHIn[2] * padfHIn[7] - padfHIn[1] * padfHIn[8]) * inv_det;
3✔
470

471
    padfHOut[0] = (padfHIn[2] * padfHIn[3] - padfHIn[0] * padfHIn[5]) * inv_det;
3✔
472
    padfHOut[3] = (padfHIn[0] * padfHIn[4] - padfHIn[1] * padfHIn[3]) * inv_det;
3✔
473
    padfHOut[6] = (padfHIn[1] * padfHIn[5] - padfHIn[2] * padfHIn[4]) * inv_det;
3✔
474

475
    return 1;
3✔
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
void *GDALCreateHomographyTransformerFromGCPs(int nGCPCount,
22✔
494
                                              const GDAL_GCP *pasGCPList)
495
{
496
    double adfHomography[9];
497

498
    if (GDALGCPsToHomography(nGCPCount, pasGCPList, adfHomography))
22✔
499
    {
500
        return GDALCreateHomographyTransformer(adfHomography);
22✔
501
    }
NEW
502
    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
void GDALDestroyHomographyTransformer(void *pTransformArg)
29✔
520

521
{
522
    if (pTransformArg == nullptr)
29✔
NEW
523
        return;
×
524

525
    HomographyTransformInfo *psInfo =
29✔
526
        static_cast<HomographyTransformInfo *>(pTransformArg);
527

528
    if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
29✔
529
    {
530
        delete psInfo;
29✔
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
int GDALHomographyTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
17,543✔
560
                            double *x, double *y, CPL_UNUSED double *z,
561
                            int *panSuccess)
562
{
563
    VALIDATE_POINTER1(pTransformArg, "GDALHomographyTransform", 0);
17,543✔
564

565
    HomographyTransformInfo *psInfo =
17,543✔
566
        static_cast<HomographyTransformInfo *>(pTransformArg);
567

568
    double *homography = bDstToSrc ? psInfo->padfReverse : psInfo->padfForward;
17,543✔
569
    int ret = TRUE;
17,543✔
570
    for (int i = 0; i < nPointCount; i++)
49,853✔
571
    {
572
        double w = homography[6] + x[i] * homography[7] + y[i] * homography[8];
32,310✔
573
        if (abs(w) < 1.0e-15)
32,310✔
574
        {
NEW
575
            panSuccess[i] = FALSE;
×
NEW
576
            ret = FALSE;
×
577
        }
578
        else
579
        {
580
            double wx =
32,310✔
581
                homography[0] + x[i] * homography[1] + y[i] * homography[2];
32,310✔
582
            double wy =
32,310✔
583
                homography[3] + x[i] * homography[4] + y[i] * homography[5];
32,310✔
584
            x[i] = wx / w;
32,310✔
585
            y[i] = wy / w;
32,310✔
586
            panSuccess[i] = TRUE;
32,310✔
587
        }
588
    }
589

590
    return ret;
17,543✔
591
}
592

593
/************************************************************************/
594
/*                 GDALSerializeHomographyTransformer()                 */
595
/************************************************************************/
596

597
CPLXMLNode *GDALSerializeHomographyTransformer(void *pTransformArg)
8✔
598

599
{
600
    VALIDATE_POINTER1(pTransformArg, "GDALSerializeHomographyTransformer",
8✔
601
                      nullptr);
602

603
    HomographyTransformInfo *psInfo =
8✔
604
        static_cast<HomographyTransformInfo *>(pTransformArg);
605

606
    CPLXMLNode *psTree =
607
        CPLCreateXMLNode(nullptr, CXT_Element, "HomographyTransformer");
8✔
608

609
    /* -------------------------------------------------------------------- */
610
    /*      Attach Homography.                                              */
611
    /* -------------------------------------------------------------------- */
612
    char szWork[300] = {};
8✔
613

614
    CPLsnprintf(
8✔
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
    CPLCreateXMLElementAndValue(psTree, "Homography", szWork);
8✔
621

622
    return psTree;
8✔
623
}
624

625
/************************************************************************/
626
/*                     GDALDeserializeHomography()                      */
627
/************************************************************************/
628

629
static void GDALDeserializeHomography(const char *pszH, double adfHomography[9])
7✔
630
{
631
    CPLsscanf(pszH, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf", adfHomography + 0,
7✔
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
void *GDALDeserializeHomographyTransformer(CPLXMLNode *psTree)
7✔
642

643
{
644
    /* -------------------------------------------------------------------- */
645
    /*        Homography                                                    */
646
    /* -------------------------------------------------------------------- */
647
    double padfForward[9];
648
    if (CPLGetXMLNode(psTree, "Homography") != nullptr)
7✔
649
    {
650
        GDALDeserializeHomography(CPLGetXMLValue(psTree, "Homography", ""),
7✔
651
                                  padfForward);
652

653
        /* -------------------------------------------------------------------- */
654
        /*      Generate transformation.                                        */
655
        /* -------------------------------------------------------------------- */
656
        void *pResult = GDALCreateHomographyTransformer(padfForward);
7✔
657

658
        return pResult;
7✔
659
    }
NEW
660
    return nullptr;
×
661
}
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc