Coverage Report

Created: 2025-06-13 06:29

/src/gdal/ogr/ogrcircularstring.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  OpenGIS Simple Features Reference Implementation
4
 * Purpose:  The OGRCircularString geometry class.
5
 * Author:   Even Rouault, even dot rouault at spatialys dot com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2010, 2014, Even Rouault <even dot rouault at spatialys dot
9
 *com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13
14
#include "cpl_port.h"
15
#include "ogr_geometry.h"
16
17
#include <cmath>
18
#include <cstring>
19
20
#include <algorithm>
21
#include <limits>
22
#include <vector>
23
24
#include "cpl_error.h"
25
#include "ogr_core.h"
26
#include "ogr_geometry.h"
27
#include "ogr_p.h"
28
29
static inline double dist(double x0, double y0, double x1, double y1)
30
0
{
31
0
    return std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0));
32
0
}
33
34
/************************************************************************/
35
/*              OGRCircularString( const OGRCircularString& )           */
36
/************************************************************************/
37
38
/**
39
 * \brief Copy constructor.
40
 *
41
 * Note: before GDAL 2.1, only the default implementation of the constructor
42
 * existed, which could be unsafe to use.
43
 *
44
 * @since GDAL 2.1
45
 */
46
47
0
OGRCircularString::OGRCircularString(const OGRCircularString &) = default;
48
49
/************************************************************************/
50
/*                  operator=( const OGRCircularString& )               */
51
/************************************************************************/
52
53
/**
54
 * \brief Assignment operator.
55
 *
56
 * Note: before GDAL 2.1, only the default implementation of the operator
57
 * existed, which could be unsafe to use.
58
 *
59
 * @since GDAL 2.1
60
 */
61
62
OGRCircularString &OGRCircularString::operator=(const OGRCircularString &other)
63
0
{
64
0
    if (this != &other)
65
0
    {
66
0
        OGRSimpleCurve::operator=(other);
67
0
    }
68
0
    return *this;
69
0
}
70
71
/************************************************************************/
72
/*                          getGeometryType()                           */
73
/************************************************************************/
74
75
OGRwkbGeometryType OGRCircularString::getGeometryType() const
76
77
12
{
78
12
    if ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED))
79
6
        return wkbCircularStringZM;
80
6
    else if (flags & OGR_G_MEASURED)
81
0
        return wkbCircularStringM;
82
6
    else if (flags & OGR_G_3D)
83
5
        return wkbCircularStringZ;
84
1
    else
85
1
        return wkbCircularString;
86
12
}
87
88
/************************************************************************/
89
/*                          getGeometryName()                           */
90
/************************************************************************/
91
92
const char *OGRCircularString::getGeometryName() const
93
94
18
{
95
18
    return "CIRCULARSTRING";
96
18
}
97
98
/************************************************************************/
99
/*                               clone()                                */
100
/************************************************************************/
101
102
OGRCircularString *OGRCircularString::clone() const
103
104
0
{
105
0
    return new (std::nothrow) OGRCircularString(*this);
106
0
}
107
108
/************************************************************************/
109
/*                           importFromWkb()                            */
110
/*                                                                      */
111
/*      Initialize from serialized stream in well known binary          */
112
/*      format.                                                         */
113
/************************************************************************/
114
115
OGRErr OGRCircularString::importFromWkb(const unsigned char *pabyData,
116
                                        size_t nSize, OGRwkbVariant eWkbVariant,
117
                                        size_t &nBytesConsumedOut)
118
119
0
{
120
0
    OGRErr eErr = OGRSimpleCurve::importFromWkb(pabyData, nSize, eWkbVariant,
121
0
                                                nBytesConsumedOut);
122
0
    if (eErr == OGRERR_NONE)
123
0
    {
124
0
        if (!IsValidFast())
125
0
        {
126
0
            empty();
127
0
            return OGRERR_CORRUPT_DATA;
128
0
        }
129
0
    }
130
0
    return eErr;
131
0
}
132
133
/************************************************************************/
134
/*                            exportToWkb()                             */
135
/*                                                                      */
136
/*      Build a well known binary representation of this object.        */
137
/************************************************************************/
138
139
OGRErr
140
OGRCircularString::exportToWkb(unsigned char *pabyData,
141
                               const OGRwkbExportOptions *psOptions) const
142
143
0
{
144
0
    if (!IsValidFast())
145
0
    {
146
0
        return OGRERR_FAILURE;
147
0
    }
148
149
0
    OGRwkbExportOptions sOptions(psOptions ? *psOptions
150
0
                                           : OGRwkbExportOptions());
151
152
    // Does not make sense for new geometries, so patch it.
153
0
    if (sOptions.eWkbVariant == wkbVariantOldOgc)
154
0
        sOptions.eWkbVariant = wkbVariantIso;
155
0
    return OGRSimpleCurve::exportToWkb(pabyData, &sOptions);
156
0
}
157
158
/************************************************************************/
159
/*                           importFromWkt()                            */
160
/*                                                                      */
161
/*      Instantiate from well known text format.  Currently this is     */
162
/*      `CIRCULARSTRING [Z] ( x y [z], x y [z], ...)',                  */
163
/************************************************************************/
164
165
OGRErr OGRCircularString::importFromWkt(const char **ppszInput)
166
167
18
{
168
18
    const OGRErr eErr = OGRSimpleCurve::importFromWkt(ppszInput);
169
18
    if (eErr == OGRERR_NONE)
170
14
    {
171
14
        if (!IsValidFast())
172
2
        {
173
2
            empty();
174
2
            return OGRERR_CORRUPT_DATA;
175
2
        }
176
14
    }
177
16
    return eErr;
178
18
}
179
180
/************************************************************************/
181
/*                            exportToWkt()                             */
182
/************************************************************************/
183
184
std::string OGRCircularString::exportToWkt(const OGRWktOptions &opts,
185
                                           OGRErr *err) const
186
0
{
187
0
    if (!IsValidFast())
188
0
    {
189
0
        if (err)
190
0
            *err = OGRERR_FAILURE;
191
0
        return std::string();
192
0
    }
193
194
0
    OGRWktOptions optsModified(opts);
195
0
    optsModified.variant = wkbVariantIso;
196
0
    return OGRSimpleCurve::exportToWkt(optsModified, err);
197
0
}
198
199
/************************************************************************/
200
/*                             get_Length()                             */
201
/*                                                                      */
202
/*      For now we return a simple euclidean 2D distance.               */
203
/************************************************************************/
204
205
double OGRCircularString::get_Length() const
206
207
0
{
208
0
    double dfLength = 0.0;
209
0
    for (int i = 0; i < nPointCount - 2; i += 2)
210
0
    {
211
0
        const double x0 = paoPoints[i].x;
212
0
        const double y0 = paoPoints[i].y;
213
0
        const double x1 = paoPoints[i + 1].x;
214
0
        const double y1 = paoPoints[i + 1].y;
215
0
        const double x2 = paoPoints[i + 2].x;
216
0
        const double y2 = paoPoints[i + 2].y;
217
0
        double R = 0.0;
218
0
        double cx = 0.0;
219
0
        double cy = 0.0;
220
0
        double alpha0 = 0.0;
221
0
        double alpha1 = 0.0;
222
0
        double alpha2 = 0.0;
223
0
        if (OGRGeometryFactory::GetCurveParameters(
224
0
                x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
225
0
        {
226
0
            dfLength += fabs(alpha2 - alpha0) * R;
227
0
        }
228
0
        else
229
0
        {
230
0
            dfLength += dist(x0, y0, x2, y2);
231
0
        }
232
0
    }
233
0
    return dfLength;
234
0
}
235
236
/************************************************************************/
237
/*                       ExtendEnvelopeWithCircular()                   */
238
/************************************************************************/
239
240
void OGRCircularString::ExtendEnvelopeWithCircular(
241
    OGREnvelope *psEnvelope) const
242
0
{
243
0
    if (!IsValidFast() || nPointCount == 0)
244
0
        return;
245
246
    // Loop through circular portions and determine if they include some
247
    // extremities of the circle.
248
0
    for (int i = 0; i < nPointCount - 2; i += 2)
249
0
    {
250
0
        const double x0 = paoPoints[i].x;
251
0
        const double y0 = paoPoints[i].y;
252
0
        const double x1 = paoPoints[i + 1].x;
253
0
        const double y1 = paoPoints[i + 1].y;
254
0
        const double x2 = paoPoints[i + 2].x;
255
0
        const double y2 = paoPoints[i + 2].y;
256
0
        double R = 0.0;
257
0
        double cx = 0.0;
258
0
        double cy = 0.0;
259
0
        double alpha0 = 0.0;
260
0
        double alpha1 = 0.0;
261
0
        double alpha2 = 0.0;
262
0
        if (OGRGeometryFactory::GetCurveParameters(
263
0
                x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
264
0
        {
265
0
            if (std::isnan(alpha0) || std::isnan(alpha2))
266
0
            {
267
0
                CPLError(CE_Failure, CPLE_AppDefined,
268
0
                         "GetCurveParameters returned NaN");
269
0
                continue;
270
0
            }
271
0
            int quadrantStart =
272
0
                static_cast<int>(std::floor(alpha0 / (M_PI / 2)));
273
0
            int quadrantEnd = static_cast<int>(std::floor(alpha2 / (M_PI / 2)));
274
0
            if (quadrantStart > quadrantEnd)
275
0
            {
276
0
                std::swap(quadrantStart, quadrantEnd);
277
0
            }
278
            // Transition trough quadrants in counter-clock wise direction.
279
0
            for (int j = quadrantStart + 1; j <= quadrantEnd; ++j)
280
0
            {
281
0
                switch ((j + 8) % 4)
282
0
                {
283
0
                    case 0:
284
0
                        psEnvelope->MaxX = std::max(psEnvelope->MaxX, cx + R);
285
0
                        break;
286
0
                    case 1:
287
0
                        psEnvelope->MaxY = std::max(psEnvelope->MaxY, cy + R);
288
0
                        break;
289
0
                    case 2:
290
0
                        psEnvelope->MinX = std::min(psEnvelope->MinX, cx - R);
291
0
                        break;
292
0
                    case 3:
293
0
                        psEnvelope->MinY = std::min(psEnvelope->MaxY, cy - R);
294
0
                        break;
295
0
                    default:
296
0
                        CPLAssert(false);
297
0
                        break;
298
0
                }
299
0
            }
300
0
        }
301
0
    }
302
0
}
303
304
/************************************************************************/
305
/*                            getEnvelope()                             */
306
/************************************************************************/
307
308
void OGRCircularString::getEnvelope(OGREnvelope *psEnvelope) const
309
310
0
{
311
0
    OGRSimpleCurve::getEnvelope(psEnvelope);
312
0
    ExtendEnvelopeWithCircular(psEnvelope);
313
0
}
314
315
/************************************************************************/
316
/*                            getEnvelope()                             */
317
/************************************************************************/
318
319
void OGRCircularString::getEnvelope(OGREnvelope3D *psEnvelope) const
320
321
0
{
322
0
    OGRSimpleCurve::getEnvelope(psEnvelope);
323
0
    ExtendEnvelopeWithCircular(psEnvelope);
324
0
}
325
326
/************************************************************************/
327
/*                     OGRCircularString::segmentize()                  */
328
/************************************************************************/
329
330
bool OGRCircularString::segmentize(double dfMaxLength)
331
0
{
332
0
    if (!IsValidFast())
333
0
        return false;
334
0
    if (nPointCount == 0)
335
0
        return true;
336
337
    // So as to make sure that the same line followed in both directions
338
    // result in the same segmentized line.
339
0
    if (paoPoints[0].x < paoPoints[nPointCount - 1].x ||
340
0
        (paoPoints[0].x == paoPoints[nPointCount - 1].x &&
341
0
         paoPoints[0].y < paoPoints[nPointCount - 1].y))
342
0
    {
343
0
        reversePoints();
344
0
        const bool bRet = segmentize(dfMaxLength);
345
0
        reversePoints();
346
0
        return bRet;
347
0
    }
348
349
0
    std::vector<OGRRawPoint> aoRawPoint;
350
0
    std::vector<double> adfZ;
351
0
    bool bRet = true;
352
0
    for (int i = 0; i < nPointCount - 2; i += 2)
353
0
    {
354
0
        const double x0 = paoPoints[i].x;
355
0
        const double y0 = paoPoints[i].y;
356
0
        const double x1 = paoPoints[i + 1].x;
357
0
        const double y1 = paoPoints[i + 1].y;
358
0
        const double x2 = paoPoints[i + 2].x;
359
0
        const double y2 = paoPoints[i + 2].y;
360
0
        double R = 0.0;
361
0
        double cx = 0.0;
362
0
        double cy = 0.0;
363
0
        double alpha0 = 0.0;
364
0
        double alpha1 = 0.0;
365
0
        double alpha2 = 0.0;
366
367
0
        aoRawPoint.emplace_back(x0, y0);
368
0
        if (padfZ)
369
0
            adfZ.emplace_back(padfZ[i]);
370
371
        // We have strong constraints on the number of intermediate points
372
        // we can add.
373
374
0
        if (OGRGeometryFactory::GetCurveParameters(
375
0
                x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
376
0
        {
377
            // It is an arc circle.
378
0
            const double dfSegmentLength1 = fabs(alpha1 - alpha0) * R;
379
0
            const double dfSegmentLength2 = fabs(alpha2 - alpha1) * R;
380
0
            if (dfSegmentLength1 > dfMaxLength ||
381
0
                dfSegmentLength2 > dfMaxLength)
382
0
            {
383
0
                const double dfVal =
384
0
                    1 + 2 * std::floor(dfSegmentLength1 / dfMaxLength / 2.0);
385
0
                if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 ||
386
0
                    std::isnan(dfVal))
387
0
                {
388
0
                    CPLError(CE_Failure, CPLE_AppDefined,
389
0
                             "segmentize nIntermediatePoints invalid: %lf",
390
0
                             dfVal);
391
0
                    bRet = false;
392
0
                    break;
393
0
                }
394
0
                const int nIntermediatePoints = static_cast<int>(dfVal);
395
0
                const double dfStep =
396
0
                    (alpha1 - alpha0) / (nIntermediatePoints + 1);
397
0
                for (int j = 1; j <= nIntermediatePoints; ++j)
398
0
                {
399
0
                    double alpha = alpha0 + dfStep * j;
400
0
                    const double x = cx + R * cos(alpha);
401
0
                    const double y = cy + R * sin(alpha);
402
0
                    aoRawPoint.emplace_back(x, y);
403
0
                    if (padfZ)
404
0
                    {
405
0
                        const double z = padfZ[i] + (padfZ[i + 1] - padfZ[i]) *
406
0
                                                        (alpha - alpha0) /
407
0
                                                        (alpha1 - alpha0);
408
0
                        adfZ.emplace_back(z);
409
0
                    }
410
0
                }
411
0
            }
412
0
            aoRawPoint.emplace_back(x1, y1);
413
0
            if (padfZ)
414
0
                adfZ.emplace_back(padfZ[i + 1]);
415
416
0
            if (dfSegmentLength1 > dfMaxLength ||
417
0
                dfSegmentLength2 > dfMaxLength)
418
0
            {
419
0
                const double dfVal =
420
0
                    1 + 2 * std::floor(dfSegmentLength2 / dfMaxLength / 2.0);
421
0
                if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 ||
422
0
                    std::isnan(dfVal))
423
0
                {
424
0
                    CPLError(CE_Failure, CPLE_AppDefined,
425
0
                             "segmentize nIntermediatePoints invalid 2: %lf",
426
0
                             dfVal);
427
0
                    bRet = false;
428
0
                    break;
429
0
                }
430
0
                int nIntermediatePoints = static_cast<int>(dfVal);
431
0
                const double dfStep =
432
0
                    (alpha2 - alpha1) / (nIntermediatePoints + 1);
433
0
                for (int j = 1; j <= nIntermediatePoints; ++j)
434
0
                {
435
0
                    const double alpha = alpha1 + dfStep * j;
436
0
                    const double x = cx + R * cos(alpha);
437
0
                    const double y = cy + R * sin(alpha);
438
0
                    aoRawPoint.emplace_back(x, y);
439
0
                    if (padfZ)
440
0
                    {
441
0
                        const double z =
442
0
                            padfZ[i + 1] + (padfZ[i + 2] - padfZ[i + 1]) *
443
0
                                               (alpha - alpha1) /
444
0
                                               (alpha2 - alpha1);
445
0
                        adfZ.emplace_back(z);
446
0
                    }
447
0
                }
448
0
            }
449
0
        }
450
0
        else
451
0
        {
452
            // It is a straight line.
453
0
            const double dfSegmentLength1 = dist(x0, y0, x1, y1);
454
0
            const double dfSegmentLength2 = dist(x1, y1, x2, y2);
455
0
            if (dfSegmentLength1 > dfMaxLength ||
456
0
                dfSegmentLength2 > dfMaxLength)
457
0
            {
458
0
                const double dfVal =
459
0
                    1 + 2 * std::ceil(dfSegmentLength1 / dfMaxLength / 2.0);
460
0
                if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 ||
461
0
                    std::isnan(dfVal))
462
0
                {
463
0
                    CPLError(CE_Failure, CPLE_AppDefined,
464
0
                             "segmentize nIntermediatePoints invalid 2: %lf",
465
0
                             dfVal);
466
0
                    bRet = false;
467
0
                    break;
468
0
                }
469
0
                int nIntermediatePoints = static_cast<int>(dfVal);
470
0
                for (int j = 1; j <= nIntermediatePoints; ++j)
471
0
                {
472
0
                    aoRawPoint.emplace_back(
473
0
                        x0 + j * (x1 - x0) / (nIntermediatePoints + 1),
474
0
                        y0 + j * (y1 - y0) / (nIntermediatePoints + 1));
475
0
                    if (padfZ)
476
0
                        adfZ.emplace_back(padfZ[i] +
477
0
                                          j * (padfZ[i + 1] - padfZ[i]) /
478
0
                                              (nIntermediatePoints + 1));
479
0
                }
480
0
            }
481
482
0
            aoRawPoint.emplace_back(x1, y1);
483
0
            if (padfZ)
484
0
                adfZ.emplace_back(padfZ[i + 1]);
485
486
0
            if (dfSegmentLength1 > dfMaxLength ||
487
0
                dfSegmentLength2 > dfMaxLength)
488
0
            {
489
0
                const double dfVal =
490
0
                    1 + 2 * std::ceil(dfSegmentLength2 / dfMaxLength / 2.0);
491
0
                if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 ||
492
0
                    std::isnan(dfVal))
493
0
                {
494
0
                    CPLError(CE_Failure, CPLE_AppDefined,
495
0
                             "segmentize nIntermediatePoints invalid 3: %lf",
496
0
                             dfVal);
497
0
                    bRet = false;
498
0
                    break;
499
0
                }
500
0
                const int nIntermediatePoints = static_cast<int>(dfVal);
501
502
0
                for (int j = 1; j <= nIntermediatePoints; ++j)
503
0
                {
504
0
                    aoRawPoint.emplace_back(
505
0
                        x1 + j * (x2 - x1) / (nIntermediatePoints + 1),
506
0
                        y1 + j * (y2 - y1) / (nIntermediatePoints + 1));
507
0
                    if (padfZ)
508
0
                        adfZ.emplace_back(padfZ[i + 1] +
509
0
                                          j * (padfZ[i + 2] - padfZ[i + 1]) /
510
0
                                              (nIntermediatePoints + 1));
511
0
                }
512
0
            }
513
0
        }
514
0
    }
515
0
    aoRawPoint.push_back(paoPoints[nPointCount - 1]);
516
0
    if (padfZ)
517
0
        adfZ.push_back(padfZ[nPointCount - 1]);
518
519
0
    CPLAssert(aoRawPoint.empty() ||
520
0
              (aoRawPoint.size() >= 3 && (aoRawPoint.size() % 2) == 1));
521
0
    if (padfZ)
522
0
    {
523
0
        CPLAssert(adfZ.size() == aoRawPoint.size());
524
0
    }
525
526
    // Is there actually something to modify?
527
0
    if (nPointCount < static_cast<int>(aoRawPoint.size()))
528
0
    {
529
0
        nPointCount = static_cast<int>(aoRawPoint.size());
530
0
        paoPoints = static_cast<OGRRawPoint *>(
531
0
            CPLRealloc(paoPoints, sizeof(OGRRawPoint) * nPointCount));
532
0
        memcpy(paoPoints, &aoRawPoint[0], sizeof(OGRRawPoint) * nPointCount);
533
0
        if (padfZ)
534
0
        {
535
0
            padfZ = static_cast<double *>(
536
0
                CPLRealloc(padfZ, sizeof(double) * aoRawPoint.size()));
537
0
            memcpy(padfZ, &adfZ[0], sizeof(double) * nPointCount);
538
0
        }
539
0
    }
540
0
    return bRet;
541
0
}
542
543
/************************************************************************/
544
/*                               Value()                                */
545
/*                                                                      */
546
/*      Get an interpolated point at some distance along the curve.     */
547
/************************************************************************/
548
549
void OGRCircularString::Value(double dfDistance, OGRPoint *poPoint) const
550
551
0
{
552
0
    if (dfDistance < 0)
553
0
    {
554
0
        StartPoint(poPoint);
555
0
        return;
556
0
    }
557
558
0
    double dfLength = 0;
559
560
0
    for (int i = 0; i < nPointCount - 2; i += 2)
561
0
    {
562
0
        const double x0 = paoPoints[i].x;
563
0
        const double y0 = paoPoints[i].y;
564
0
        const double x1 = paoPoints[i + 1].x;
565
0
        const double y1 = paoPoints[i + 1].y;
566
0
        const double x2 = paoPoints[i + 2].x;
567
0
        const double y2 = paoPoints[i + 2].y;
568
0
        double R = 0.0;
569
0
        double cx = 0.0;
570
0
        double cy = 0.0;
571
0
        double alpha0 = 0.0;
572
0
        double alpha1 = 0.0;
573
0
        double alpha2 = 0.0;
574
575
        // We have strong constraints on the number of intermediate points
576
        // we can add.
577
578
0
        if (OGRGeometryFactory::GetCurveParameters(
579
0
                x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
580
0
        {
581
            // It is an arc circle.
582
0
            const double dfSegLength = fabs(alpha2 - alpha0) * R;
583
0
            if (dfSegLength > 0)
584
0
            {
585
0
                if ((dfLength <= dfDistance) &&
586
0
                    ((dfLength + dfSegLength) >= dfDistance))
587
0
                {
588
0
                    const double dfRatio =
589
0
                        (dfDistance - dfLength) / dfSegLength;
590
591
0
                    const double alpha =
592
0
                        alpha0 * (1 - dfRatio) + alpha2 * dfRatio;
593
0
                    const double x = cx + R * cos(alpha);
594
0
                    const double y = cy + R * sin(alpha);
595
596
0
                    poPoint->setX(x);
597
0
                    poPoint->setY(y);
598
599
0
                    if (getCoordinateDimension() == 3)
600
0
                        poPoint->setZ(padfZ[i] * (1 - dfRatio) +
601
0
                                      padfZ[i + 2] * dfRatio);
602
603
0
                    return;
604
0
                }
605
606
0
                dfLength += dfSegLength;
607
0
            }
608
0
        }
609
0
        else
610
0
        {
611
            // It is a straight line.
612
0
            const double dfSegLength = dist(x0, y0, x2, y2);
613
0
            if (dfSegLength > 0)
614
0
            {
615
0
                if ((dfLength <= dfDistance) &&
616
0
                    ((dfLength + dfSegLength) >= dfDistance))
617
0
                {
618
0
                    const double dfRatio =
619
0
                        (dfDistance - dfLength) / dfSegLength;
620
621
0
                    poPoint->setX(paoPoints[i].x * (1 - dfRatio) +
622
0
                                  paoPoints[i + 2].x * dfRatio);
623
0
                    poPoint->setY(paoPoints[i].y * (1 - dfRatio) +
624
0
                                  paoPoints[i + 2].y * dfRatio);
625
626
0
                    if (getCoordinateDimension() == 3)
627
0
                        poPoint->setZ(padfZ[i] * (1 - dfRatio) +
628
0
                                      padfZ[i + 2] * dfRatio);
629
630
0
                    return;
631
0
                }
632
633
0
                dfLength += dfSegLength;
634
0
            }
635
0
        }
636
0
    }
637
638
0
    EndPoint(poPoint);
639
0
}
640
641
/************************************************************************/
642
/*                          CurveToLine()                               */
643
/************************************************************************/
644
645
OGRLineString *
646
OGRCircularString::CurveToLine(double dfMaxAngleStepSizeDegrees,
647
                               const char *const *papszOptions) const
648
0
{
649
0
    OGRLineString *poLine = new OGRLineString();
650
0
    poLine->assignSpatialReference(getSpatialReference());
651
652
0
    const bool bHasZ = getCoordinateDimension() == 3;
653
0
    for (int i = 0; i < nPointCount - 2; i += 2)
654
0
    {
655
0
        OGRLineString *poArc = OGRGeometryFactory::curveToLineString(
656
0
            paoPoints[i].x, paoPoints[i].y, padfZ ? padfZ[i] : 0.0,
657
0
            paoPoints[i + 1].x, paoPoints[i + 1].y, padfZ ? padfZ[i + 1] : 0.0,
658
0
            paoPoints[i + 2].x, paoPoints[i + 2].y, padfZ ? padfZ[i + 2] : 0.0,
659
0
            bHasZ, dfMaxAngleStepSizeDegrees, papszOptions);
660
0
        poLine->addSubLineString(poArc, (i == 0) ? 0 : 1);
661
0
        delete poArc;
662
0
    }
663
0
    return poLine;
664
0
}
665
666
/************************************************************************/
667
/*                        IsValidFast()                                 */
668
/************************************************************************/
669
670
OGRBoolean OGRCircularString::IsValidFast() const
671
672
14
{
673
14
    if (nPointCount == 1 || nPointCount == 2 ||
674
14
        (nPointCount >= 3 && (nPointCount % 2) == 0))
675
2
    {
676
2
        CPLError(CE_Failure, CPLE_NotSupported,
677
2
                 "Bad number of points in circular string : %d", nPointCount);
678
2
        return FALSE;
679
2
    }
680
12
    return TRUE;
681
14
}
682
683
/************************************************************************/
684
/*                            IsValid()                                 */
685
/************************************************************************/
686
687
OGRBoolean OGRCircularString::IsValid() const
688
689
0
{
690
0
    return IsValidFast() && OGRGeometry::IsValid();
691
0
}
692
693
/************************************************************************/
694
/*                         hasCurveGeometry()                           */
695
/************************************************************************/
696
697
OGRBoolean
698
OGRCircularString::hasCurveGeometry(int /* bLookForNonLinear */) const
699
12
{
700
12
    return TRUE;
701
12
}
702
703
/************************************************************************/
704
/*                         getLinearGeometry()                        */
705
/************************************************************************/
706
707
OGRGeometry *
708
OGRCircularString::getLinearGeometry(double dfMaxAngleStepSizeDegrees,
709
                                     const char *const *papszOptions) const
710
0
{
711
0
    return CurveToLine(dfMaxAngleStepSizeDegrees, papszOptions);
712
0
}
713
714
//! @cond Doxygen_Suppress
715
/************************************************************************/
716
/*                     GetCasterToLineString()                          */
717
/************************************************************************/
718
719
static OGRLineString *CasterToLineString(OGRCurve *poGeom)
720
0
{
721
0
    CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
722
0
             poGeom->getGeometryName());
723
0
    delete poGeom;
724
0
    return nullptr;
725
0
}
726
727
OGRCurveCasterToLineString OGRCircularString::GetCasterToLineString() const
728
0
{
729
0
    return ::CasterToLineString;
730
0
}
731
732
/************************************************************************/
733
/*                        GetCasterToLinearRing()                       */
734
/************************************************************************/
735
736
static OGRLinearRing *CasterToLinearRing(OGRCurve *poGeom)
737
0
{
738
0
    CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
739
0
             poGeom->getGeometryName());
740
0
    delete poGeom;
741
0
    return nullptr;
742
0
}
743
744
OGRCurveCasterToLinearRing OGRCircularString::GetCasterToLinearRing() const
745
0
{
746
0
    return ::CasterToLinearRing;
747
0
}
748
749
//! @endcond
750
751
/************************************************************************/
752
/*                            IsFullCircle()                            */
753
/************************************************************************/
754
755
int OGRCircularString::IsFullCircle(double &cx, double &cy,
756
                                    double &square_R) const
757
0
{
758
0
    if (getNumPoints() == 3 && get_IsClosed())
759
0
    {
760
0
        const double x0 = getX(0);
761
0
        const double y0 = getY(0);
762
0
        const double x1 = getX(1);
763
0
        const double y1 = getY(1);
764
0
        cx = (x0 + x1) / 2;
765
0
        cy = (y0 + y1) / 2;
766
0
        square_R = (x1 - cx) * (x1 - cx) + (y1 - cy) * (y1 - cy);
767
0
        return TRUE;
768
0
    }
769
    // Full circle defined by 2 arcs?
770
0
    else if (getNumPoints() == 5 && get_IsClosed())
771
0
    {
772
0
        double R_1 = 0.0;
773
0
        double cx_1 = 0.0;
774
0
        double cy_1 = 0.0;
775
0
        double alpha0_1 = 0.0;
776
0
        double alpha1_1 = 0.0;
777
0
        double alpha2_1 = 0.0;
778
0
        double R_2 = 0.0;
779
0
        double cx_2 = 0.0;
780
0
        double cy_2 = 0.0;
781
0
        double alpha0_2 = 0.0;
782
0
        double alpha1_2 = 0.0;
783
0
        double alpha2_2 = 0.0;
784
0
        if (OGRGeometryFactory::GetCurveParameters(
785
0
                getX(0), getY(0), getX(1), getY(1), getX(2), getY(2), R_1, cx_1,
786
0
                cy_1, alpha0_1, alpha1_1, alpha2_1) &&
787
0
            OGRGeometryFactory::GetCurveParameters(
788
0
                getX(2), getY(2), getX(3), getY(3), getX(4), getY(4), R_2, cx_2,
789
0
                cy_2, alpha0_2, alpha1_2, alpha2_2) &&
790
0
            fabs(R_1 - R_2) < 1e-10 && fabs(cx_1 - cx_2) < 1e-10 &&
791
0
            fabs(cy_1 - cy_2) < 1e-10 &&
792
0
            (alpha2_1 - alpha0_1) * (alpha2_2 - alpha0_2) > 0)
793
0
        {
794
0
            cx = cx_1;
795
0
            cy = cy_1;
796
0
            square_R = R_1 * R_1;
797
0
            return TRUE;
798
0
        }
799
0
    }
800
0
    return FALSE;
801
0
}
802
803
/************************************************************************/
804
/*                       get_AreaOfCurveSegments()                      */
805
/************************************************************************/
806
807
//! @cond Doxygen_Suppress
808
double OGRCircularString::get_AreaOfCurveSegments() const
809
0
{
810
0
    double dfArea = 0.0;
811
0
    for (int i = 0; i < getNumPoints() - 2; i += 2)
812
0
    {
813
0
        const double x0 = getX(i);
814
0
        const double y0 = getY(i);
815
0
        const double x1 = getX(i + 1);
816
0
        const double y1 = getY(i + 1);
817
0
        const double x2 = getX(i + 2);
818
0
        const double y2 = getY(i + 2);
819
0
        double R = 0.0;
820
0
        double cx = 0.0;
821
0
        double cy = 0.0;
822
0
        double alpha0 = 0.0;
823
0
        double alpha1 = 0.0;
824
0
        double alpha2 = 0.0;
825
0
        if (OGRGeometryFactory::GetCurveParameters(
826
0
                x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
827
0
        {
828
            // Should be <= PI in absolute value.
829
0
            const double delta_alpha01 = alpha1 - alpha0;
830
0
            const double delta_alpha12 = alpha2 - alpha1;  // Same.
831
            // http://en.wikipedia.org/wiki/Circular_segment
832
0
            dfArea += 0.5 * R * R *
833
0
                      fabs(delta_alpha01 - sin(delta_alpha01) + delta_alpha12 -
834
0
                           sin(delta_alpha12));
835
0
        }
836
0
    }
837
0
    return dfArea;
838
0
}
839
840
//! @endcond
841
842
/************************************************************************/
843
/*                           get_Area()                                 */
844
/************************************************************************/
845
846
double OGRCircularString::get_Area() const
847
0
{
848
0
    if (IsEmpty() || !get_IsClosed())
849
0
        return 0;
850
851
0
    double cx = 0.0;
852
0
    double cy = 0.0;
853
0
    double square_R = 0.0;
854
855
0
    if (IsFullCircle(cx, cy, square_R))
856
0
    {
857
0
        return M_PI * square_R;
858
0
    }
859
860
    // Optimization for convex rings.
861
0
    if (IsConvex())
862
0
    {
863
        // Compute area of shape without the circular segments.
864
0
        double dfArea = get_LinearArea();
865
866
        // Add the area of the spherical segments.
867
0
        dfArea += get_AreaOfCurveSegments();
868
869
0
        return dfArea;
870
0
    }
871
872
0
    OGRLineString *poLS = CurveToLine();
873
0
    const double dfArea = poLS->get_Area();
874
0
    delete poLS;
875
876
0
    return dfArea;
877
0
}
878
879
/************************************************************************/
880
/*                        get_GeodesicArea()                            */
881
/************************************************************************/
882
883
double OGRCircularString::get_GeodesicArea(
884
    const OGRSpatialReference *poSRSOverride) const
885
0
{
886
0
    if (IsEmpty())
887
0
        return 0;
888
889
0
    if (!get_IsClosed())
890
0
    {
891
0
        CPLError(CE_Failure, CPLE_AppDefined, "Non-closed geometry");
892
0
        return -1;
893
0
    }
894
895
0
    if (!poSRSOverride)
896
0
        poSRSOverride = getSpatialReference();
897
898
0
    auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
899
0
    return poLS->get_GeodesicArea(poSRSOverride);
900
0
}
901
902
/************************************************************************/
903
/*                        get_GeodesicLength()                          */
904
/************************************************************************/
905
906
double OGRCircularString::get_GeodesicLength(
907
    const OGRSpatialReference *poSRSOverride) const
908
0
{
909
0
    if (IsEmpty())
910
0
        return 0;
911
912
0
    if (!poSRSOverride)
913
0
        poSRSOverride = getSpatialReference();
914
915
0
    auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
916
0
    return poLS->get_GeodesicLength(poSRSOverride);
917
0
}
918
919
//! @cond Doxygen_Suppress
920
921
/************************************************************************/
922
/*                           ContainsPoint()                            */
923
/************************************************************************/
924
925
int OGRCircularString::ContainsPoint(const OGRPoint *p) const
926
0
{
927
0
    double cx = 0.0;
928
0
    double cy = 0.0;
929
0
    double square_R = 0.0;
930
0
    if (IsFullCircle(cx, cy, square_R))
931
0
    {
932
0
        const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
933
0
                                   (p->getY() - cy) * (p->getY() - cy);
934
0
        return square_dist < square_R;
935
0
    }
936
0
    return -1;
937
0
}
938
939
/************************************************************************/
940
/*                       IntersectsPoint()                              */
941
/************************************************************************/
942
943
int OGRCircularString::IntersectsPoint(const OGRPoint *p) const
944
0
{
945
0
    double cx = 0.0;
946
0
    double cy = 0.0;
947
0
    double square_R = 0.0;
948
0
    if (IsFullCircle(cx, cy, square_R))
949
0
    {
950
0
        const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
951
0
                                   (p->getY() - cy) * (p->getY() - cy);
952
0
        return square_dist <= square_R;
953
0
    }
954
0
    return -1;
955
0
}
956
957
//! @endcond